CC BY
34
УДК 550.388.2
ВОССТАНОВЛЕНИЕ ТРЕХМЕРНОГО РАСПРЕДЕЛЕНИЯ ЭЛЕКТРОННОЙ КОНЦЕНТРАЦИИ ИОНОСФЕРЫ ПО ИОНОГРАММАМ DYNASONDE-21
© 2009 г. Г.А. Жбанков, Н.А. Заботин
Научно-исследовательский институт физики Южного федерального университета, 344090, г. Ростов-на-Дону, пр. Стачки, 194, zhbankov@ip.rsu.ru
Research Institute of Physics of Southern Federal University, 344090, Rostov-on-Don, Stachki Ave, 194, zhbankov@ip.rsu.ru
Задача восстановления реального высотного профиля электронной концентрации (инверсии) по данным ионозонда имеет историю, начинающуюся с Эпплтона. Она всегда формулировалось как прикладная математическая задача и поэтому требовала набора существенных физических идеализаций. Главное упрощение, свойственное всем методам восстановления профиля — предположение о плоской горизонтально-стратифицированной ионосфере, сводящее проблему только к одному измерению.
Ключевые слова: ионосфера, распространение радиоволн, неоднородности, электронная концентрация, ионограмма.
The problem of electron density inversion of digital ionograms is reconsidered from the viewpoints of new possibilities and of modern requirements. The data processing system of an advanced ionosonde (the dynasonde) provides accurate measurements not only of echo group range but also of direction of arrival, among other physical parameters, thus yielding the three-dimensional distribution of apparent echoloca-tions in each ionogram recording. An iterative ray-tracing approach is described here to recover the parameters of a quite sophisticated three-dimensional (so-called wedge-stratified ionosphere) model of the local electron density distribution, character!/ing its actual vertical Ne(h} profile together with horizontal gradients and general tilts. This approach is implemented in the algorithm introduced here, is named "NeXtYZ".
Keywords: ionosphere, wave propagation, disturbance, electron density, ionograms.
Свойства реальной ионосферы значительно отличаются от идеализированной картины [1]. Фактически все процессы в ионосфере создают горизонтальные и наклонные градиенты электронной концентрации. К ним можно отнести солнечный терминатор, акустические гравитационные волны, перемещающиеся ионосферные возмущения, ионосферные бури, экваториальную аномалию, плазменные пузыри, выпадающие частицы и общую широтно-долготную зависимость ионосферы. В комбинации эти эффекты заставляют отраженный сигнал распространяться по неканоническим траекториям, которые только случайно являются «вертикальными», т.е. фактические данные находятся
обычно в более или менее явном разногласии с обычными предположениями.
Современные цифровые ионосферные зонды предоставляют возможность точного измерения и времени группового распространения, и направления прихода для каждого отраженного луча на ионограмме. Явные трехмерные представления результатов радиозондирования представлены на рис. 1.
Каждая точка на этих графиках представляет конец «вектора группового пути» [2], построенного в соответствии с направлением прихода луча, измеренным диназондом, эквивалентна длине половины времени группового распространения, которое определено с большой точностью методом стационарной фазы.
УС кт
XL hm
б
Рис. 1. Пример трехмерного представления результатов радиозондирования, полученных автоматизированной процедурой анализа данных измерений dynasonde-21 DSND в EISCAT Tromso (а) и в обсерватории Bear Lake, штат Юта (б)
52
200,0
100,0
а
Фактические положения точек отражения (где выполняются условия отражения) не измерены непосредственно и будут отличаться из-за рефракции в неоднородной магнитоактивной ионосферной плазме. В «прямой задаче» эти эффекты могут быть рассчитаны по лучевой траектории, если известно трехмерное распределение ионосферной плазмы. В процессе восстановления профиля расчет лучевой траектории становится компонентом в процедуре как звено, связывающее пространственное распределение электронной плотности с измеренными значениями отраженных лучей и их угловыми координатами. Это - основная стратегия новой схемы инверсии NeXtYZ.
Задача трехмерной инверсии может быть сформулирована как восстановление параметров параметрической модели, которая описывает в локальном масштабе и вертикальные, и горизонтальные градиенты плотности ионосферной плазмы. Модель клиновидно-стратифицированной ионосферы (КСИ) - соответствующая замена для прежней модели плоско-стратифицированной ионосферы.
Клиновидно-стратифицированная ионосфера
В модели КСИ поверхности постоянной плотности электронной концентрации задаются последовательно с небольшими приращениями по плазменной частоте /р{ в диапазоне высот Н{ по вертикальной оси; наклон каждой структурной плоскости характеризуется двумя горизонтальными компонентами вектора нормали к ней пх{, пу{. Нормаль к поверхности электронной плотности определяет локальное направление общего градиента в слое.
Модель КСИ предполагает «клиновидное» распределение электронной плотности между двумя соседними структурными плоскостями. Основное требование заключается в том, что плазменная частота является постоянной на каждой такой плоскости. Поверхности структуры не должны пересекаться в пределах «поля зрения» ионозонда. Так как КСИ - локальная модель распределения плотности электронной концентрации, пересечения структурных плоскостей допустимы вне области ее применимости.
Стоит заметить, что распространение зондирующих сигналов в модели КСИ является существенно трехмерным. Аналитические приближения этого процесса не могут дать качественного приближения, так что при расчете траекторий необходимо использовать численные методы.
Основной особенностью модели КСИ является возможность описания изменения трехмерного наклона локальных ионосферных слоев с высотой. Опция одномерного (т.е. вертикального) профиля не потеряна и является частным случаем.
Версии алгоритма, действующие в настоящее время, воспроизводят параметры индивидуального набора плоскостей КСИ для каждой ионограммы. Такой подход применим для многих целей, что дает возможность для развития, включая более тщательное описание сложных ионосферных ситуаций.
Основные принципы алгоритма инверсии
Рассмотрим предварительные замечания о влиянии ионосферных неоднородностей на параметры зондирующих сигналов, а также некоторые требования, предъявляемые этим к искомому алгоритму инверсии электронной плотности.
1. Статистические особенности и неопределенность данных: алгоритм должен легко обрабатывать данные, содержащие существенный случайный компонент или неопределенность, вызванную естественной пространственной изменчивостью в ионосферной электронной плотности ниже выбранного уровня расчета.
2. Всесторонний ввод: алгоритм должен использовать всю доступную экспериментальную информацию о пространственных свойствах распределения электронного плотности.
3. Физически значимое пространственное разрешение: конкретные ионосферные условия должны определять форму профиля электронной плотности на текущем этапе расчета; подход, основанный на приспособлении заданной комбинации функций ко всему ионосферному профилю, неприемлем.
Чтобы удовлетворять этим требованиям, NeXtYZ реализует следующие принципы и шаги.
1. Большинство упрощающих предположений современных алгоритмов инверсии заменено численным расчетом лучевых траекторий для всех отобранных лучей на ионограмме, приняты во внимание их наблюдаемые углы прихода.
2. Инверсия начинается с начального уровня, продолжаясь последовательностью тонких слоев плазменной частоты («клиньев»).
3. Инверсия проводится на основе метода наименьших квадратов, минимизируя разницу между расчетными и фактическими временами распространения и направлениями для индивидуального отражения.
4. Имея дело с неизбежной «начальной проблемой» (практические наблюдения никогда не начинаются с наименьших значений плазменной частоты), первый шаг решается одновременно для параметров самого низкого плазменного слоя (клина) и для трех поправочных коэффициентов модели ионизации Тит-териджа [3-5]. Эта модель сшивается с начальным интервалом инверсии.
5. Все плазменные клинья, относящиеся к наблюдаемой монотонной части ионосферного слоя, рассматриваются, используя одну и ту же повторяющуюся процедуру. Параметры всех нижележащих клиньев остаются неизменными. Вертикальное положение йг+1 следующей структурной плоскости (соответствующей плазменной частоте /рМ) и ее
ориентация (nx
yi+\
) определяются в результате
специальной процедуры минимизации для разности между измеренными и расчетными параметрами того набора лучей, которые отражаются между
и /Рм (Рис- 2).
6. В случае, если наблюдения дают признаки Е-Б долины, ее общая форма определяется из ранее использованной модели [6]. Параметры модели долины
подбираются в соответствии с лучами, имеющими отношение к F-области по той же самой процедуре оптимизации, которая используется при получении характеристик следующего плазменного клина.
7. В области пиков Е- и F-слоев проблема минимизации решается для параметров модели слоя Чепмена вместо расчета вертикальных положений независимых структурных плоскостей. Ориентации плоскостей также определяются по специальной процедуре оптимизации.
8. Расчетные положения точек отражения и времени распространения для каждого луча, использованного в процедуре инверсии, сохраняются при выходе из процедуры.
Рис. 2. Геометрическое представление метода оптимизации. Каждая итерация процедуры требует нахождения групповых времен распространения для всех лучей, отраженных в пределах текущего клина. Лучевые траектории и соответствующие им точки отражения схематически представлены линиями в двумерной проекции. Также показаны два ранее определенных клина и текущий искомый клин (1+1) со своей верхней границей
Численный расчет лучевых траекторий
Алгоритм численного расчета лучевых траекторий, используемый NeXtYZ, базируется в настоящее время на методе характеристик [7] для уравнения эйконала вида
1
Н&г,/,® у^ \
2-kW
(1)
Основываясь на этом, получим систему дифференциальных уравнений, описывающих поведение луча, волнового вектора и времени распространения в виде
йг _ ж _р_1дп2
йт др 2 ф
йг дг 2 дг '
й 2 + 1 а 2 йт 2 да
где г - независимая вспомогательная переменная; г = ^.у.г у _ радиус-вектор в декартовой сис-_к |к„
волновой вектор; |к0| = к0 = ю/с ; п - показатель
преломления падающей волны в магнитоактивной плазме (заданный классической формулой Эпплтона); t - время распространения; со - радиочастота. Для случая распространения сигнала в Е- и F-областях столкновениями в плазме можно пренебречь.
Процедура оптимизации методом наименьших квадратов
Клиновидно-стратифицированная модель ионосферы требует только стандартного и относительно простого распределения электронной плотности в пределах единичных клиньев. Предположим, что параметры всех основных клиньев были уже найдены, т.е. существуют три параметра, которые полностью управляют этим распределением для текущего клина: вертикальное положение /г;+1 его верхней границы (следующей структурной плоскости, соответствующей плазменной частоте и ориентация этой
плоскости (характеризуемая горизонтальными компонентами вектора ее нормали пхМ и пуп]). Множество лучей, отраженных в пределах каждого клина, обычно находится в соответствии с процедурами анализа данных диназонда. Для целей инверсии мы требуем наличия, по крайней мере, 10 таких лучей. Если это необходимо, стандартный шаг по плазменной частоте, используемый алгоритмом NeXtYZ, для некоторых клиньев может увеличиваться.
Каждый отраженный сигнал характеризуется набором физических параметров, определенных в соответствии с предварительной процедурой анализа данных диназонда. К ним относятся амплитуда; поляризация; два угла прихода в; и ср 1 : групповой путь
R
0 j
измеренный методом стационарной фазы; ус-
теме координат; р = -.—г = Д - нормированный
редненная фаза Я * ^; оценки соответствующих им
ошибок. Все параметры, исключая фазу, используются NeXtYZ. Значение амплитуды в децибелах применяется как весовой коэффициент при некоторых усредняющих операциях и процедурах выборки. Поляризация («обыкновенная» или «необыкновенная») определяет соответствующее выражение для показателя преломления падающей волны в магнитоактив-ной плазме при траекторных расчетах луча. Углы прихода обеспечивают начальные условия на поверхности для системы уравнений луча и определяют две других полезных характеристики, средние углы прихода (0;+1 и Ф;+1) для лучей, отраженных в пределах клина. Значение группового пути используется без каких-либо изменений в процедуре оптимизации.
Точки отражения лучей в пределах клина могут распределяться произвольно. Если распределение электронной плотности достаточно гладкое, точки отражения имеют тенденцию располагаться вдоль линии, не обязательно прямой, независимой переменной для которой может быть частота или время или их комбинация. В присутствии мелкомасштабных структур (очень типичная ситуация) точки отражения рассеяны в пространстве (как и показано на рис. 1). NeXtYZ не налагает искусственных ограничений на
пространственное распределение точек отражения, при условии, что траектория луча рассчитана достаточно точно. Если рассчитанное положение точки отражения луча оказывается ниже текущего клина, то этот луч исключается из процесса оптимизации.
И групповой путь , и углы прихода 0) и
луча, отраженного в пределах текущего клина, чувствительны к параметрам текущего клина (йг+1, пхМ и п м). Цель и главная идея процедуры оптимизации
состоят в том, чтобы определить наилучшие значения //, |, пхп | и пуМ , минимизируя (последовательно, а
не совместно) две величины:
1) отклонение Мг'+1 = б'+и 3 для Рас"
считанного (р) ) и измеренного (Л'.) значений групповых времен;
2) расстояние между расположением зонда (точкой отсчета системы координат) и «точкой возврата» для луча, излученного из этой точки под углами ©г+1 и Фг+1, определенными выше.
Второй компонент процедуры оптимизации увеличивает ее чувствительность к наклону ионосферных слоев; это оправдано в соответствии с требованием обеспечения условия наклонной структуры ионосферы, при котором лучи возвращаются в точку излучения. В реальной ионосфере мелкомасштабные неоднородности облегчают согласование с этим условием для всего набора обнаруженных отражений. Существующая версия клиновидно-стратифицированной модели ионосферы не включает никаких мелкомасштабных структур, и автоматическое возвращение к точке расположения зонда не гарантируется для всех отражений, так что мы применяем условие минимизации к лучу со средним значением направления прихода для лучей в пределах клина.
Как отмечено, два компонента процедуры оптимизации явно не объединяются в NeXtYZ; мы считаем более эффективным их поочередное осуществление в последовательных шагах итерационного процесса для текущего клина. При таком условии достигается высокий уровень стабильности обеих частей процедуры оптимизации, что делает ненужным применение различных специальных методов регуляризации.
Оценки погрешности инверсии
Набор основных параметры отражений, передаваемые диназондом в процедуру минимизации, избыточен. Это дает возможность оценить погрешности представления реального ионосферного слоя результирующей КСИ моделью.
Погрешность оценивает неопределенность результата работы самой процедуры инверсии метода NeXtYZ. Качество оптимизации определяется остаточным значением минимизируемой величины
Л1<!', | = 4?', 11 у - Л,' | у ^ , относящейся к / -му клину. Чем меньше , тем решение более точное. Однако минимальное значение остатка фактически никогда не достигает нуля. Это можно рассматривать
как статистическую характеристику остаточного локального (1 -го клина) несоответствия между набором данных и решением. Это - готовая оценка локальной погрешности инверсии. Необходимо только преобразовать эту величину из виртуального значения в реальную высоту.
Проверка и тестирование NeXtYZ в режиме реального времени
Полная версия NeXtYZ-алгоритма, описанного в предыдущих разделах, была проверена с модельными наборами данных. Процесс моделирования начинается с определения невозмущенного вертикального профиля (модель, состоящая из параболических сегментов). Профиль определяет вертикальные положения для последовательности структурных плоскостей в модели КСИ. Далее задается модель возмущения, приводящая к изменению положения структурных плоскостей (например, волнообразная). Неоднородности накладываются на модельное крупномасштабное распределение электронной концентрации для определенной цели - соответствовать обычному ионосферному «спокойному состоянию».
Полученные смоделированные параметры отражений передаются в NeXtYZ как входная информация. Проводя процедуру инверсии с модельными данными, мы можем сравнить ее результаты с модельным вертикальным / ^ ^ профилем и модельной высотной зависимостью двух компонентов вектора нормали к плоскости п. Тесты с использованием такой процедуры подтвердили хорошую работу №Х1У2-алгоритма. Вертикальный профиль всегда восстанавли-
вался удовлетворительно, ограничения в дневных условиях были связаны (как и для всех методов инверсии) с его зависимостью от качества модели Е-Р долины. Существенные ошибки в компонентах вектора нормали п проявляются, главным образом, там, где общий градиент электронной плотности является малым (например, около максимума Р-слоя) и, также как и у высотного профиля, в области долины, где информация о наклоне отсутствует.
Алгоритм NeXtYZ был экстенсивно проверен в режиме реального времени в работе с данными dynasonde-21, начиная с декабря 2004 г. Все другие программы обработки данных dynasonde-21 теперь базируются на результатах NeXtYZ-инверсии. В течение периода испытания результаты анализа были доступны на web-сайте http://www.ngdc.noaa.gov/stp/ЮNO/Dynasonde/. В течение длительного периода тестирования параллельные вычисления проводились параллельно методом POLAN, и результаты, полученные обоими алгоритмами, были представлены для сравнения в режиме реального времени. Эти результаты доступны для ретроспективного сравнения в базе данных dynasonde на этом же сайте.
В течение нескольких лет приблизительно 225 000 ионограмм были обработаны с помощью данного метода. Эта долгосрочная серия тестов демонстрирует другую важную особенность NeXtYZ алгоритма. Он является намного более устойчивым к грубым ошибкам выбора следов и плохим данным. Ионосферные
профили, полученные NeXtYZ, значительно более однородны и устойчивы к плохим или недостаточным данным, чем те же самые, полученные с помощью процедуры POLAN.
Оптимизация кода NeXtYZ проводилась параллельно с его развитием. К настоящему времени анализ типичных ионограмм dynasonde-21 NeXtYZ алгоритмом инверсии профиля занимает приблизительно 2,5 мин на стандартных PC с процессором Pentium-4 3,4 ГГц. Для упрощенной версии NeXtYZ-Lite среднее время анализа - всего 0,5 мин. Эта характеристика делает обе версии NeXtYZ подходящими не только для научных исследований, но также и для современного эксплуатационного применения.
Из изложенного выше следует, что горизонтальные градиенты электронной плотности присутствуют в ионосфере в широком диапазоне пространственных масштабов и амплитуд. Они имеют важное практическое значение для основных характеристик ионосферного состояния, вызывая радиомерцания, ошибки и отказы GPS, ошибки в навигационных и локационных измерениях. Они также лежат в основе изучения ионосферы и как причина, и как результат нестабильности и неоднородностей плазмы. Ионосферное зондирование с его специальной «селективностью уровня отражения» во всем частотном диапазоне является потенциально одной из самых чувствительных частей диагностики этих градиентов; NeXtYZ - первая реализация такого потенциала.
Алгоритм NeXtYZ заменяет POLAN для целей восстановления профиля электронной концентрации в прикладном аспекте. Основное развитие ядра NeXtYZ -полный трехмерный алгоритм - закончено. Он был проверен на распределениях электронной плотности, описанных КСИ и модельных ионограммах. Параметры КСИ в модельных расчетах воспроизведены удов-
летворительно. Это также было проверено с декабря 2004 г. в режиме реального времени по данным работающего устройства dynasonde-21.
Созданию NeXtYZ-модели способствовали два главных обстоятельства. Во-первых, программа анализа данных диназонда DSND обеспечивает высокую плотность, обычно несколько тысяч отражений на ионограмме, каждое снабжено точными физическими параметрами. Это позволяет решить проблему оптимизации для каждого «клина» электронной плотности с высокой точностью. Во-вторых, современные PC быстро производят численные траекторные расчеты. Обработка данных ионограмм в режиме реального времени с использованием NeXtYZ происходит в разумное время (0,5-2,5 мин) и не требует дорогих решений аппаратных средств ЭВМ.
Литература
1. Appleton E. V. Some notes on wireless methods of investigating the electrical structure of the upper atmosphere // Proc. Phys. Soc. London, 1930. P. 321.
2. PaulA.K., Wright J.W., FedorL.S. The interpretation of ionospheric drift measurements. Angle of arrival and group path (echolocation) measurements from digitized ionospheric soundings: The group path vector // J. Atmos. Terr. Phys. 1974. № 2. P. 193.
3. Titheridge J.E. Modeling the peak of the ionospheric E-layer // J. Atmos. Sol. Terr. Phys. 2000. № 1. P. 93.
4. Titheridge J.E. Production of the low-latitude night E layer // J. Geophys. Res. 2001. № 6. P. 781.
5. Titheridge J.E. Ionization below the night F2 layer: A global model // J. Atmos. Sol. Terr. Phys. 2003. № 10. P. 1035.
6. Titheridge J.E. Model results for the daytime ionospheric E and valley regions // J. Atmos. Sol. Terr. Phys. 2003. № 2. P. 129.
7. Кравцов Ю.А., Орлов Ю.И. Геометрическая оптика неоднородных сред. М., 1980. 304 с.
Поступила в редакцию_12 августа 2008 г.
