CC BY
83
УДК 550.388
РАДИОТОМОГРАФИЯ ИОНОСФЕРЫ С ПРИМЕНЕНИЕМ ВЫСОКООРБИТАЛЬНЫХ НАВИГАЦИОННЫХ СИСТЕМ
В. Е. Куницын, Е. С. Андреева, М. А. Кожарин, И. А. Нестеров
(.кафедра физики атмосферы) E-mail: kunitsyn@phys.msu.su
Предложены варианты радиотомографии ионосферы с применением высокоорбитальных навигационных систем типа GPS, ГЛОНАСС. Рассмотрение высокоорбитальной радиотомографии проведено в сопоставлении с радиотомографией на низкоорбитальных навигационных системах. Приведены результаты радиотомографических реконструкций по различным экспериментальным данным и сравнение результатов высокоорбитальной и низкоорбитальной радиотомографии.
Введение
С начала 90-х гг. интенсивно развиваются методы радиотомографии (РТ), позволяющие восстанавливать пространственную структуру распределения электронной концентрации при радиозондировании ионосферы [1, 2]. Существующие низкоорбитальные (НО) спутниковые навигационные системы (типа российских спутников «Цикада» и американских «Транзит» — Navy Navigation Satellite System), имеющие практически круговую орбиту на высоте около 1000-1150 км, и наземные цепочки приемников дают возможность получать серии томографических данных по различным лучам. В РТ-экспериментах прием двух когерентных спутниковых сигналов на частотах 150 и 400 МГц и регистрация разности фаз между ними (приведенной фазы) осуществляется на сети из нескольких наземных приемных станций, расположенных вдоль траектории спутников и на расстояниях порядка сотен километров. Такие системы позволяют получать двумерные РТ-сечения ионосферы на расстояниях в 1-3 тысячи километров за времена порядка 10-15 мин. В настоящее время вышло несколько сотен работ, связанных с РТ ионосферы, с основными результатами РТ-исследований можно ознакомиться в следующих монографиях и обзорах [1-5].
Для реализации метода РТ на базе НО навигационных систем необходимо создание сетей (линеек) приемников. В настоящее время существует около десятка таких сетей, которые активно используются для исследовательских целей. С другой стороны, в последние годы активно развиваются сети приемников высокоорбитальных (ВО) навигационных систем, таких как GPS и ГЛОНАСС. В настоящее время функционирует более тысячи приемников ВО систем. В данной статье рассмотрен вариант РТ ионосферы: РТ на ВО системах. Рассмотрение ВО РТ целесообразно провести в сопоставлении с НО РТ.
1. РТ на низкоорбитальных спутниках — НО РТ
Измерения приведенной фазы <р в нескольких приемных точках являются данными для томографической реконструкции. Интеграл от электронной концентрации N вдоль лучей между приемником на Земле и передатчиком на спутнике пропорционален абсолютной (полной) фазе Ф, включая неизвестную начальную фазу щ:
аХге У Ыйа = Ф = <р + <р0, (1)
ь
где А — длина зондирующей волны, ге — классический радиус электрона, безразмерный коэффициент пропорциональности а ~ 1 определяется выбором частот.
Перепишем уравнение (1) в операторном виде, добавив типичный некоррелированный экспериментальный шум £:
РЛ: Ф + £. (2)
здесь Р — проекционный оператор, переводящий двумерное распределение N в набор одномерных функций (проекций) Ф. Таким образом, задача томографической реконструкции состоит в решении системы линейных интегральных уравнений (2) и нахождении электронной концентрации N. Однако не существует способов непосредственного решения системы линейных интегральных уравнений вида (2). Одним из возможных вариантов является дискретизация (аппроксимация) проекционного оператора Р и замена его на дискретный оператор Ь. В этом случае получаем соответствующую систему линейных уравнений (СЛУ):
ЬЛ: Ф+ £ + /<;. Е = ЪЫ ^ РЖ, (3)
где Е — ошибка аппроксимации, коррелированная ошибка, зависящая от самого решения. При этом оператор Ь будет определяться выбором типа аппроксимации или представления искомой функ-
ции N. Отметим, что уравнения (2) и (3) эквивалентны. Но в РТ-экеперименте ошибка аппроксимации неизвестна, и мы фактически решаем другую систему:
ЬЛ; Ф + £. (4)
Система (4) неэквивалентна системе (3), и все ошибки и искажения, связанные с ошибкой аппроксимации Е, будут проявляться в реконструкции. Ошибка аппроксимации Е существенно варьируется в зависимости от типа аппроксимации. Чем ближе аппроксимация проекционного оператора Ь к истинному проекционному оператору Р, тем меньше погрешность Е. Однако никакие методы решения СЛУ не позволяют выявлять структуры, вклад которых в данные будет существенно меньше присутствующих при обработке искажений, связанных с ошибкой аппроксимации Е [6].
Важно отметить, что для решения томографической СЛУ (4) необходимо знать абсолютную (полную) фазу Ф, т.е. определить неизвестную начальную фазу (ра. Погрешности в оценке значений начальных фаз (ра у разных приемников могут привести к противоречивости и несогласованности данных и, как результат, к низкому качеству реконструкций. Поэтому был разработан метод фазоразностной РТ или РТ по разности линейных интегралов на близких лучах, не требующий определения начальной фазы (5) [1, 7, 8]. СЛУ для фазоразностной РТ (5) определяется соответствующей разностью:
ЪЫ = Ф 1
> А /V = Ф Ф' = П-£, (5)
Ъ'ы = Ф' ^
где ЬЖ = Ф — исходная СЛУ и Ь'ЛГ = Ф' — СЛУ, сформированная по набору близких лучей. Фазоразностная РТ дает заметно лучшие результаты и обладает более высокой чувствительностью по сравнению с фазовыми методами, что подтверждается экспериментальными результатами [9-12]. Разрешение задачи лучевой РТ в линейной постановке составляет 20-30 км по горизонтали и 30-40 км по вертикали. Если учитывать рефракцию [13] зондирующих лучей, то пространственное разрешение метода лучевой РТ ионосферы можно улучшить до 10-20 км.
Существует большое число алгоритмов решения систем линейных уравнений (4), (5) как прямых, так и итерационных. В настоящее время в задачах лучевой РТ ионосферы наиболее часто используются итерационные алгоритмы, также применяют и неитерационные алгоритмы: метод сингулярного разложения и его модификации, регуляризации среднеквадратичного отклонения, ортогональной декомпозиции, максимальной энтропии, квадратичного программирования и его разновидности, байесовский подход и т.д. [4, 14]. За последние пятнадцать лет в процессе моделирования и проведения большого числа экспериментальных
РТ-реконструкций были найдены эффективные комбинации различных методов и алгоритмов, дающие наиболее качественные реконструкции. Важнейшим априорным условием задачи лучевой РТ ионосферы является неотрицательность решения. Причем в итерационных алгоритмах это условие легко учесть в процессе восстановления, в отличие от неитерационных, в частности, в байесовском подходе в реконструкциях могут появляться отрицательные значения электронной концентрации. Еще одним из преимуществ итерационных алгоритмов является возможность учитывать в процессе реконструкции краевые эффекты и априорную информацию о возможной структуре ионосферы [15].
Применение различных методов РТ, типов аппроксимаций проекционного оператора и алгоритмов решения СЛУ вместе с дополнительными условиями (геометрией, выбором сетки дискретизации и т. д.) порождает большое число решений задачи лучевой РТ ионосферы (4)-(5), причем все эти решения эквивалентны экспериментальным данным с заданной точностью. Поэтому был предложен метод построения ансамбля решений [6] и нахождения среднего решения по ансамблю, который позволяет оценить погрешность решения (или «дисперсию» решения) и тем самым определелить вероятность появления тех или иных деталей реконструируемой структуры, что особенно важно для практических приложений задач лучевой РТ. При увеличении числа членов ансамбля мы получим в пределе единственное решение СЛУ (4), (5). В реальном РТ-эксперименте число членов ансамбля достаточно ограничить единицами сотен членов, и добавление новых в ансамбль в этом случае практически не будет влиять на среднее решение по ансамблю. Причем среднее решение по ансамблю является наилучшим (оптимальным) решением томографической задачи, так как в процессе усреднения различного вида искажения и артефакты компенсируются. Таким образом, построение ансамбля решений позволяет получить устойчивое решение задачи лучевой РТ.
В последние годы лучевая РТ стала одним из основных инструментов исследования распределений электронной плотности в ионосфере. Были проведены и в настоящее время реализуются серии РТ-эке-периментов в различных регионах мира: в Европе, Америке и Азии. Впервые в мире экспериментальные РТ-реконструкции были получены в марте-апреле 1990 г. сотрудниками МГУ и ПГИ РАН [9]. На рис. 1 в координатах географическая широта-высота в километрах приведено в изолиниях (в единицах 106 электрон ем-3) одно из первых полученных РТ-сечений ионосферы между Москвой и Мурманском, на котором хорошо виден провал в районе 63° северной широты. Проведенные многочисленные эксперименты показали сложность и разнообразие форм провала, причем его ширина, наклон и глубина варьируются в широких пределах.
Ь, км
700 ■
500 '
300 -
66
Широта
Рис. 1. РТ-реконструкция главного ионосферного провала по трассе Мурманск-Москва по данным от 8.04.90, 00:43 1Т
В ходе РТ-экепериментов в области приэкваториальных широт был выявлен ряд структурных особенностей экваториальной аномалии (ЭА): ориентация сформировавшегося ядра ЭА в полуденные часы вдоль направления магнитного поля Земли, существенная асимметрия краев ЭА, а также характерное сжатие ионосферы севернее ядра ЭА. Пример восстановления структуры Р- и Е-областей в районе ЭА в Юго-Восточной Азии (между Манилой и Шанхаем) приведен на рис. 2. На реконструкции хорошо видно «продавливание» нижней ионосферы плазменным потоком в области широт 26-28°. Исследования ЭА были проведены в ходе совместных работ МГУ, университетов Иллинойса, Тайваня и Уханя [16-19]. Отметим, что восстановление структуры Е-области ионосферы методами РТ существенно сложнее, поскольку вклад Е-области в регистрируемые данные значительно меньше вклада .Р-области [20]. Метод НО РТ позволяет не только получать двумерные сечения концентрации плазмы, но и определять потоки плазмы, рассматривая последовательные во времени сечения [19]. На рис. 3 дан пример определения потоков плазмы (произведения плотности и скорости плазмы), хорошо иллюстрирующий явление фонтан-эффекта.
И, км
700 ■
500
20 22
24 26 28
Широта
Рис. 2. РТ-реконструкция Е- и Е-областей северного гребня экваториальной аномалии между Шанхаем и Манилой по данным от 7.10.94, 15:40 ЦТ
к, 700
600
500
400
300
КМ
14.09.94 (12:50-13:30ЬТ)
Величина потока (109 см 2с')
\ 1 1 < '
МИ''
V I 4 / ' 1 V V
12 14 16 18 20 22 24 26 28
Широта
Рис. 3. Пример потоков ионосферной плазмы по НО
РТ-данным от 14.09.94 (12:50-13:30 1Л"), иллюстрирующий явление фонтан-эффекта
РТ-методы позволили исследовать как известные ионосферные структуры, так и найти ряд новых (квазиволновые структуры, структура ядра экваториальной аномалии, наклонный провал ионизации, пальцеобразные структуры, локализованные неоднородности и др.). Исследованы особенности структуры и динамики экваториальной аномалии и области главного ионосферного провала. Разработан метод определения потоков ионосферной плазмы по томографическим данным [19]. РТ-реконструкции перемещающихся ионосферных возмущений дали информацию о параметрах возмущений и позволили исследовать атмоеферно-ионоеферные взаимодействия. Проведены РТ-исследования сильных возмущений ионосферы, вызванных антропогенными факторами, в частности возмущений, вызванных стартами ракет [4], промышленными взрывами [21], мощным КВ излучением [4]. Методами статистической РТ были получены распределения интенсивности флук-туаций ионосферной плазмы [22].
В ходе ряда экспериментов проведены сопоставления РТ-результатов с данными ионозондов и радаров некогерентного рассеяния. Российско-американский эксперимент показал высокое качество РТ-рекон-струкций и совпадение результатов РТ и радарных сечений в рамках точностей обоих методов [11]. Аналогично следует констатировать совпадение данных ионозондов и РТ в рамках ограничений обоих методов [18].
РТ ионосферы, основанная на данных НО навигационных спутников, позволяет исследовать двумерные сечения ионосферы, а при наличии нескольких линеек приемников, расположенных на расстояниях порядка нескольких сотен километров друг от друга, и трехмерную структуру ионосферы. Пожалуй основным существенным ограничением метода НО РТ является необходимость создания систем со многими линейками приемников. С другой
стороны, в мире уже существуют сети приемников ВО навигационных систем. Отсюда естественно возникает задача РТ на базе ВО навигационных систем.
2. Пространственно-временная томография ионосферы по данным сети ВО систем
С развертыванием глобальных навигационных систем GPS и ГЛОНАСС появился новый инструмент, который позволяет непрерывно проводить измерения характеристик проходящих через ионосферу радиосигналов. Информация, постоянно получаемая сетью приемников GPS, охватывающей почти весь земной шар, дает возможность постановки задач региональной и глобальной реконструкции распределения электронной плотности в ионосфере. Основной особенностью таких задач, относящихся к типу томографических с неполными данными, является их высокая размерность. Сравнительно небольшая угловая скорость движения ВО спутников GPS приводит к необходимости учета временной изменчивости ионосферы, т.е. к постановке задачи четырехмерной (4D) томографии (три пространственные координаты и время). Особенностью 4D томографической задачи является существенная неполнота данных — лучи спутник-приемник проходят не через все точки пространства, а в областях с малым числом приемников образуются зоны отсутствия данных, что требует разработки специальных подходов для преодоления проблем, связанных с неединственностью решения. Рассмотрим два возможных подхода.
Для решения задачи 4D РТ по данным ВО систем можно использовать подход, разработанный в двумерной спутниковой ионосферной томографии. Основная идея метода заключается в представлении искомого поля плотности электронов в виде разложения по некоторым локальным базисным функциям. В этом случае набор линейных интегралов или разностей линейных интегралов преобразуется в систему линейных уравнений. Однако в отличие от двумерной НО РТ здесь необходима дополнительная процедура интерполяции найденных решений в области отсутствия данных. Реализация этого подхода в регионах с густой сетью станций (Северная Америка, Европа) дала хорошие результаты при выборе достаточно грубой сетки и подходящих сплайнов различной гладкости.
Второй подход основывается на идее поиска изначально достаточно гладких решений задачи, таких, чтобы используемые алгоритмы обеспечивали хорошую интерполяцию в области отсутствия данных. Ищется решение, минимизирующее некоторую соболевскую норму, на бесконечном множестве решений исходной (недоопределенной) томографической задачи.
2.1. Данные для высожоорбитальной радиотомографии (ВО РТ)
GPS/ГЛОНАСС приемники позволяют получать информацию различного типа. Для целей мониторинга ионосферы используются значения фазы радиосигнала (фазового пути) при прохождении им пути от спутника до наземной приемной станции для двух рабочих частот. В обеих системах эти частоты фиксированы, например для GPS эти частоты равны: Д = 1575.42 МГц и /2 = 1227.60 МГц. Соответствующие данные принято называть L\ и L2; они представляют собой фазовые пути радиосигналов, измеренные в длинах зондирующих волн. Кроме того, можно использовать псевдодальности Pi и Р2 (групповые пути радиосигналов), измеренные по времени распространения цуга волн на частотах Д и Д. Фазовые данные L\ и L2 позволяют вычислить полное электронное содержание (ТЕС — total electron content) вдоль распространения луча по следующей формуле:
TEC
. д д
/i/f
Л2-/! к
— + const,
(6)
где К = 40.308 м3 с-2, с = 3 • 108 м/с — скорость света в вакууме.
Подчеркнем, что, используя только данные фазовых измерений Ь\ и Ь2, можно вычислить ТЕС только с точностью до неизвестной константы, которая и указана в формуле в виде дополнительного аддитивного слагаемого.
Псевдодальности Р\ и Р2 также позволяют вычислить интеграл электронной концентрации по формуле
Р2-Р1
ТЕС =
К
(А 7?)
Однако такие данные гораздо более зашумлены, чем фазовые. Уровень шума составляет, как правило, 30-50% и в некоторых случаях достигает даже 100%. Кроме того, как показали наши исследования, ТЕС, вычисленный таким образом, также содержит некоторую аддитивную константу, зависящую от станции и спутника и вызванную, по-видимому, чае-тотно-завиеимыми задержками в аппаратуре. В то же время зашумленность фазовых данных составляет 0.1-0.2 TECU, т.е. менее 1%, хотя иногда может достигать нескольких процентов. Здесь TECU — единица измерения TEC (1016 м-2). Эти факты проиллюстрированы записями типичных экспериментальных данных на рис. 4.
Таким образом, использование для радиотомографии ионосферы фазовых данных представляется более предпочтительным. Для устранения неизвестного сдвига в фазовых данных (проблема «initial phase») применяются фазоразноетные методы, разработанные ранее для РТ на НО спутниках. Традиционно томографические задачи используют данные по набору линейных интегралов. Фазоразноетные мето-
(7)
100 80 60 40 20 о -20
140 120 100 80 60 40 20 О -20 -40
1050
Г, мин
О ш
2
■чО
"о
О
ш
500 550 '/. МИН
600 650 700
Т, мин
600 650 700
Т, мин
Рис. 4. Примеры зависимости экспериментальных данных — псевдодальности (кривые /) от зенитного угла спутника (сплошные линии 2) для разных приемников (а). Примеры ТЕС, вычисленные по групповым (кривые 1, 2, 3) и фазовым
(сплошные линии 4) данным для разных приемников (б)
ды позволяют решать задачи РТ, используя только разности линейных интегралов. Иными словами, в качестве входных данных для задачи РТ используются не абсолютные значения ТЕС, а производные по времени ёТЕС/Л.
2.2. Подход с В-сплайн аппроксимацией
Формулировка задачи РТ на ВО спутниках аналогична формулировке задачи для РТ на НО спутниках. Традиционно здесь исходные экспериментальные данные выражают через ТЕС (6), (7). Как показано выше, для задачи РТ целесообразно использовать ТЕС, определяемый по фазовым данным (6). Используемый нами фазоразноетный подход позволяет в качестве данных брать разности линейных интегралов или производные по времени сГГЕС/Л.
Обозначим концентрацию электронов через Ж(г, , тогда выражение для ТЕС будет иметь вид
УлГ(г,*)<й = ТЕС^, (8)
Ьг
где интеграл вычисляется вдоль г-го луча распространения радиосигнала от спутника к приемнику.
Поскольку частота сигнала порядка одного гигагерца, то искривлением луча при прохождении им через ионосферу и атмосферу можно пренебречь и считать, что Ь^ — это прямая линия. Таким образом, для решения поставленной задачи требуется найти решение приведенной выше системы линейных интегральных уравнений. В данной работе мы использовали подход, разработанный в двумерной спутниковой ионосферной томографии. Основная идея метода заключается в представлении искомого поля плотности электронов в виде разложения по некоторым базисным функциям, т.е.
ЛГ(г,*) = £с^(г,*). (9)
з
В нашей работе в качестве базисных функций использованы В -сплайны степени с нулевой по третью. Записывая интеграл электронной концентрации вдоль луча и вычисляя производную по времени, получим:
и э и
Решая эту систему относительно коэффициентов C'j, получим приближение к искомой функции N(r,t). Для решения системы применяются итерационные методы, используемые в лучевой радиотомографии ионосферы: ART, DART, MART, SIRT и др. [4, 14].
На практике при получении РТ-реконструкций использовались комбинации алгоритмов. За последние пятнадцать лет в процессе моделирования и валидации различных вариантов НО РТ были найдены лучшие для данных условий (экспериментальных и состояния ионосферы) комбинации методов алгоритмов. Частично этот богатый опыт НО РТ может быть использован и в случае ВО РТ. Однако существенные отличия метода ВО РТ делают необходимым проведение серий исследований и поиска эффективных методов и алгоритмов. В случае трехмерной ВО РТ в отличие от двумерной РТ на НО спутниках необходима дополнительная процедура интерполяции найденных решений в области отсутствия данных. Здесь была также использована сплайн-интерполяция решений.
2.3. Подход со сглаживанием
Второй подход основывается на идее поиска достаточно гладких решений задачи, таких, чтобы используемые алгоритмы обеспечивали хорошую интерполяцию в областях отсутствия данных. Например, выберем некоторую соболевскую норму и будем искать решение, минимизирующее эту норму на бесконечном множестве решений исходной (недо-определенной) томографической задачи AN = D:
mm ||/-/о||^.
На пути практической реализации этого подхода встречаются сложности, связанные с решением минимизационной задачи с ограничениями. Прямой путь ее решения, основанный на использовании метода неопределенных множителей Лагранжа, приводит к СЛУ с матрицами высокой размерности (ввиду большого числа лучей), не обладающими специальной структурой, которая могла бы облегчить решение СЛУ. Поэтому мы предложили иной, итерационный, метод решения той же минимизационной задачи, являющийся вариацией метода S1RT с дополнительным сглаживаем путем фильтрации итерационных приращений по пространственным переменным. Метод позволяет учитывать априорную информацию, выраженную как в начальном приближении для итераций, так и в виде весовых коэффициентов, задающий относительные масштабы вариаций электронной концентрации на различных высотах.
Как известно, если система линейных уравнений Ах = у является недоопределенной, то итерационные методы типа ART
Уг ~ (aj,xfc)
xfc+i=xfc
(Щ, ai)
-a i
(здесь аj — строки матрицы А, т. е. (а¿)j = Ау) или S1RT
v с®*'®*)
(где А¡¡,pi — параметры релаксации) сходятся к нормальному решению относительно начального приближения хо, т. е. к решению задачи
min ||х-хо||2,
Ах=у
где норма является евклидовой.
Модифицируем итерационные методы так, чтобы они сходились к нормальному решению относительно иной нормы. Заметим, что минимизируемая норма порождается скалярным произведением, используемым в формулах, задающих итерации. Используем скалярное произведение, соответствующее некоторой матрице 5:
(z,x)s = (5z, 5х).
Потребуем, чтобы матрица 5*5 была обратимой. Тогда
(z,x) = ((S*S)-4x).
Перепишем итерационные формулы в терминах нового скалярного произведения:
xfc+i=xfc
Vi ~ Hxfc)sa, (ai'aDs
(iM).
-a,
где а^ = (5*5) ^-зц. Такие итерационные процессы будут сходиться к решению с минимальной нормой ЦхЦз = ^(ж, 5*5х), что и требовалось получить. В частности, ЭШТ с весами
(а
Pi =
з
можно записать в виде
х*+1 х/.' + ^Б)-1 £а, (у* - (а,, х*)).
г
Параметр г) для ускорения сходимости можно выбирать так, чтобы невязка была минимальной.
Вычислительная эффективность алгоритма зависит от простоты решения системы уравнений с матрицей 5*5. Для соболевских норм эта матрица будет иметь ленточную структуру. В случае равномерных прямоугольных сеток для решения этого уравнения может быть применен метод Фурье с использованием быстрого преобразования Фурье (БПФ), что требует порядка М 1о§2 М операций, здесь М — число узлов сетки.
Важной особенностью задач мониторинга ионосферы и магнитосферы при использовании ВО спутников является сферичность земной поверхности
с радиусом R и ионосферных слоев, что обусловливает использование сферической системы координат и соответствующих ей сеток для метода конечных элементов. Эффективность изложенной выше схемы решения задачи с существенно неполными данными основана на использовании БПФ для обращения разностного аналога дифференциального оператора, возникающего при минимизации производных искомой функции. При использовании БПФ для обращения необходимо, чтобы дифференциальный оператор был «однородным», т.е. его форма не зависела от выбора системы координат. В сферической системе координат оператор Лапласа таковым не является. При рассмотрении пространственно-временной постановки задачи, например
mini J ^2 + A2(v/)2+r2/^ dvdt^
в выражение для вариационной производной минимизируемого функционала добавляется еще и слагаемое, пропорциональное второй производной по времени:
/ - (frr2)r + r2cos2eUv +
rz COS в J
при этом слагаемые / и r2ftt от координат не зависят, А, г — характерные масштабы пространства и времени соответственно.
Такой «неоднородный» оператор, содержащий коэффициенты, зависящие от координат, можно эффективно обращать в одномерном случае, когда матрица соответствующего разностного оператора имеет п-диагональную структуру (в случае производной второго порядка матрица будет трехдиагональной и задача решается методом прогонки). Поэтому достаточно привести уравнение к виду, в котором все коэффициенты зависели бы только от одной какой-либо координаты, и тогда станет возможным использование следующего алгоритма: вначале по «однородным» координатам (г, <р) производим преобразование Фурье, заменяя производные степенями новых координат, затем обращаем оператор по оставшейся координате, после чего совершаем обратное преобразование Фурье.
В качестве «неоднородной» координаты оставим переменную 0, зависимость от которой наиболее сложная, а от радиальной избавимся путем перехода к переменой р = In и небольшого изменения исходной постановки задачи, введя в минимизируемый функционал весовую функцию, зависящую от радиуса:
mini J j> + A2^(v/)2 + r2/^ dvdt
Теперь вариационная производная минимизируемого функционала равна
А2 / 1
/ - д2 [fpp + 3/р + ^¡20^ +
а уравнение для нахождения / после преобразования Фурье по переменным р, <р, т примет вид
^"Бг ~~Ц (¡есжв) =д.
Rz cosy V /в
Введение для пространственных производных весовой функции r2/R2 приведет к более сильному сглаживанию решения в областях, удаленных от поверхности Земли, что можно считать вполне оправданным, так как вариации электронной концентрации в магнитосфере существенно меньше, чем в ионосфере (из-за того, что меньше сама электронная концентрация).
Получившиеся после дискретизации системы уравнений решаются методом прогонки. На основе приведенной схемы был реализован алгоритм решения задачи томографии по данным GPS. Как моделирование процесса томографической реконструкции, так и реконструкции при использовании реальных данных показали успешную работу алгоритма.
3. Моделирование и примеры
томографических реконструкций
Для оценки качества РТ-реконструкций на базе ВО систем было проведено компьютерное моделирование. Реконструировались возможные распределения электронной концентрации, различные неоднородности, структуры типа уединенных волн. Результаты моделирования показали приемлемое качество восстановления для квазистатических структур, мало меняющихся на протяжении одного часа. Хотя разрешение ВО РТ значительно хуже разрешения НО РТ, как правило, горизонтальное разрешение не лучше 100 км в Европе и на большей части территории США. Только в районе Южной Калифорнии возможно достичь разрешения в 30-50 км. Вертикальное разрешение при этом также около 100 км. Разрешение определяется плотностью и параметрами приемной сети. На рис. 5, а показана европейская сеть GPS приемников и траектории движения (рис. 5, б) подионоеферных точек (точек пересечения высоты 300 км лучом епутник-прием-ник) за один час. Как уже отмечалось выше, лучи спутник-приемник проходят не через все точки пространства, и в областях с малым числом приемников образуются обширные области отсутствия данных.
При моделировании применялись различные комбинации подходов и алгоритмов. Учитывая многолетний опыт НО РТ, можно отметить, что для по-
Рис. 5. Европейская сеть GPS приемников (а); траектории движения подионосферных точек за один час (б)
иска эффективных методов и алгоритмов требуются достаточно длительные исследования. В ближайшем будущем предстоит установить эффективные методы предобработки сигналов и лучшие алгоритмы для заданных классов решений.
Для иллюстрации результатов разработанных методов приведем пример реконструкции ионосферы над Европой для периода сильнейшей магнитной бури (Кр = 9) 28-31 октября 2003 г. [23] путем 4D томографической реконструкции данных GPS на европейской сети приемников в интервале [—10°, 45°]
по долготе и [30°, 80°] по широте с часовой дискретизацией по времени. По данным 41) реконструкции на рис. 6, 7 в изолиниях приведены результаты расчета электронной концентрации в единицах 106 ем-3 на высоте 300 км во время сильных возмущений в 23:00 иТ 30 октября и в 02:00 ОТ 31 октября 2003 г. Как уже сообщалось ранее [23], над исследуемым регионом располагалось пятно высокой ионизации. Здесь приведено положение пятна на входе и выходе из европейского региона. Полученные значения электронной концентрации в центре ноч-
Рис. 6. Реконструкция распределения электронной концентрации на высоте 300 км 30.10.2003 в 23:00 UT
по ВО РТ-данным
Рис. 7. Реконструкция распределения электронной концентрации на высоте 300 км 31.10.2003 в 02:00 UT
по ВО РТ-данным
ного пятна (2-2.5) • 106 ем-3 заметно превышают типичные дневные значения. Вертикальные сечения электронной концентрации ионосферы вдоль нулевого меридиана в 00:00 UT и 01:00 UT 31 октября 2003 г. показаны на рис. 8, 9. Область повышенной концентрации, расположенная между широтами 50° и 55° N, отражает вертикальную структуру упомянутого пятна [23]. Конечно, разрешающая способность и качество данных ВО РТ реконструкций заметно уступают реконструкциям НО РТ, однако основные, существенные для геофизики структуры, типа провала или пятна восстанавливаются неплохо. Ранее было проведено сопоставление критических частот слоя F по данным 4D РТ с данными ряда европейских ионозондов. Сопоставление результатов реконструкции с независимыми измерениями ионозондов показало [23] хорошие качество и эффективность разработанных подходов и алгоритмов 4D РТ по данным GPS.
На рис. 10 приведен пример ВО РТ реконструкции другой важной ионосферной структуры — экваториальной аномалии. Здесь результаты 4D реконструкции представлены распределением полного электронного содержания в единицах TECU. На реконструкции хорошо видны северный и южный гребни экваториальной аномалии [16-19].
На верхних половинах рис. 11, 12 представлены меридиональные НО РТ сечения ионосферы по данным от 21 октября 2002 г., 18:34 UT и 22 октября 2002 г., 12:38 UT в районе Аляски. Приемники
Широта
Рис. 8. Реконструкция вертикального сечения электронной концентрации ионосферы вдоль нулевого меридиана 31.10.2003 в 00:00 ЦТ по ВО РТ-данным
Широта
Рис. 9. Реконструкция вертикального сечения электронной концентрации ионосферы вдоль нулевого меридиана 31.10.2003 в 01:00 ЦТ по ВО РТ-данным
22.02.2004 5:00 UT
110 120 130 140 150 160
Долгота
Рис. 10. Реконструкция ТЕС в области экваториальной аномалии 22.02.2004, 05:00 UT по ВО РТ-данным
h, км
800
70 72
Широта
Рис. 11. Меридиональные сечения ионосферы в районе Аляски, полученные с помощью НО РТ по данным от 22.10.2002, 18:34 ЦТ (а) и ВО РТ по данным от 22.10.2002, 18:30 ЦТ (б)
располагались в Cordova (60.495° N, 15.17 Е), Gakona (62.399° N, ^145.157° Е) и Delta (63.902° N, — Mo. 2-1 Е). Экспериментальные РТ-данные доступны на сайте [24]. На нижних половинах рис. 11, 12 представлены меридиональные сечения, полученные с помощью ВО РТ. Видно, что квазирегулярная ионосфера (рис. 11) восстанавливается с высокой точностью. Ионосфера с возмущениями восстанавливается существенно хуже, ВО РТ явно не до-
h, км
800
56 58 60 62 64 66 68 70
68 70
Широта
Рис. 12. Меридиональные сечения ионосферы в районе Аляски, полученные с помощью НО РТ по данным от 22.10.2002, 12:38 ЦТ (а) и ВО РТ по данным от 22.10.2002, 12:30 ЦТ (б)
стает разрешения и не все существующие максимумы проявляются на реконструкции. Отметим, что на НО РТ-реконструкции возмущения типа АГВ (рис. 12, а) наблюдаются за 23 ч до землетрясения с эпицентром (63.514° N. — М 7.912 Е), которое произошло 23 октября 2002 г. в 11:27 иТ (магнитуда М = 6.36, глубина порядка 4.2 км). Данные примеры относятся к магнитоспокойному периоду, когда Кр варьировалось в интервале от 1 до 2.
Заключение
НО РТ позволяет реконструировать «мгновенные» (5-8 мин) 2D сечения электронной концентрации над заданным регионом. Интервал между реконструкциями зависит от числа действующих спутников и составляет в настоящее время 30-100 мин. Тем самым можно реконструировать 2D сечения не-однородностей, волновых и квазиволновых структур (в частности, ионосферные следы АГВ), структур типа уединенных волн.
ВО РТ позволяет реконструировать 4D распределения электронной концентрации (3D реконструкции каждый час или 0.5 ч). Это дает возможность реконструировать 3D распределения электронной плотности в регионе, восстанавливать структуру крупных неоднородностей, включая движущиеся неоднородности типа уединенных волн. ВО РТ дает возможность определять энергетику волновых возмущений в регионе, а также проводить региональный мониторинг уровня электронной концентрации. Однако пространственная разрешающая способность ВО РТ заметно ниже, чем в случае НО РТ, так же как и качество реконструкций волновых и квазиволновых структур.
НО РТ и ВО РТ прекрасно дополняют друг друга для целей мониторинга ионосферы. Комбинация двух систем позволит реализовать эффективный 4D региональный мониторинг ионосферы.
Авторы благодарны Северо-Западной Исследовательской Ассоциации (NWRA) за экспериментальные данные по Аляске [24].
Работа выполнена при финансовой поддержке РФФИ (гранты 05-05-65145, 04-05-64671) и 1NTAS (грант 03-50-4872).
Литература
1. Куницын В.Е., Терещенко Е.Д. Томография ионосферы. М., 1991.
2. Kunitsyn V.E., Tereshchenko E.D. // IEEE Antennas Propagation Magn. 1992. 34 P. 22.
3. Leitinger R. 11 Rev. Radio Sei. 1999. P. 581.
4. Kunitsyn V.E., Tereshchenko E.D. Ionospheric Tomography. Springer-Verlag, 2003.
5. Pryse S.E. /У Surveys in Geophysics. 2003. 24. P. 1.
6. Andreeisa E. S., Franke J.S., Kunitsyn V.E. et al. // Radio Sei. 2001. 36. P. 299.
7. Kunitsyn V.E., Andreeva E.S., Razinkov O.G., Tereshchenko E.D. // Int. J. Imaging Syst. Technol. 1994. 5. P. 128.
8. Андреева E.C., Куницын B.E. Терещенко Е.Д. // Геомагни-тизм и аэрономия. 1992. 10. С. 849.
9. Андреева Е.С., Куницын В.Е., Терещенко Е.Д. и др. // Письма в ЖЭТФ. 1990. 52. С. 145.
10. Kunitsyn V.E., Andreeva E.S., Tereshchenko E.D. et al. // Int. J. Imaging Syst. Technol. 1994. 5. P. 112.
11. Foster J.C., Kunitsyn V.E., Tereshenko E.D. et al. // Ibid. P. 148.
12. Kunitsyn V.E., Tereshchenko E.D., Andreeva E.S. et al. // Ann. Geophys. 1995. 13. P. 1242.
13. Андреева E.C., Куницын B.E., Попов А.Ю. // Вестн. Моск. ун-та. Физ. Астрон. 1999. №6. С. 21 (Moscow University Phys. Bull. 1999. N 6. P. 24).
14. Kunitsyn V.E., Andreeva E.S., Popov A.Yu., Razinkov O.G. 11 Ann. Geophys. 1995. 13. P. 1263.
15. Андреева E.C. // Радиотехн. и электроника. 2004. 49, № 1. С. 5.
16. Andreeva Е. S., Franke J.S., Yeh К.С., Kunitsyn V.E. // Geophys. Res. Lett. 2000. 27. P. 2465.
17. Yeh К.С., Franke S.J., Andreeva E.S., Kunitsyn V.E. // Ibid. 2001. 28. P. 4517.
18. Franke S.J., Yeh K.C., Andreeva E.S., Kunitsyn V.E. // Radio Sei. 2003. 38. P. 11-1.
19. Kunitsyn V.E., Andreeva E.S., Franke S.J., Yeh K.C. // Geophys. Res. Lett. 2003. 30. P. 1851.
20. Андреева E.C. // Вестн. Моск. ун-та. Физ. Астрон. 2004. №2. С. 62 (Moscow University Phys. Bull. 2004. N 2. P. ?).
21. Андреева E.C., Гохберг М.Б., Куницын B.E. и др. // Кос-мич. исследования. 2001. 39, №1. Р. 13.
22. Tereshchenko Е. D., Kozlova М.О., Kunitsyn V.E., Andreeva. E.S. 11 Radio Sei. 2004. 39. P. RS1S35-1.
23. Куницын В.E., Кожарин М.А., Нестеров И.А., Козлова М.О. И Вестн. Моск. ун-та. Физ. Астрон. 2004. №6. С. 67 (Moscow University Phys. Bull. 2004. N 6).
24. http://www.haarp.alaska.edu/haarp/data.fcgi.
Поступила в редакцию 09.12.04
