1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 63, № 5, 2023

Журнал вычислительной математики и математической физики, 2023, T. 63, № 5, стр. 803-826

Усвоение данных для двумерного уравнения амбиполярной диффузии в модели ионосферы земли

В. П. Дымников 1, 2, Д. В. Кулямин 1, 2, П. А. Останин 1, 3, *, В. П. Шутяев 1, 3

1 Институт вычислительной математики им. Г.И. Марчука РАН
119333 Москва, ул. Губкина, 8, Россия

2 Институт прикладной геофизики им. Е.К. Федорова
129128 Москва, ул. Ростокинская, 9, Россия

3 Московский физико-технический институт (национальный исследовательский университет)
141701 М.о., Долгопрудный, Институтский пер., 9, Россия

* E-mail: ostanin.pavel@phystech.edu

Поступила в редакцию 11.08.2022
После доработки 05.10.2022
Принята к публикации 02.02.2023

Полный текст (PDF)

Аннотация

В работе рассматривается задача вариационного усвоения данных наблюдений для двумерной диффузионной модели F слоя ионосферы Земли ИВМ РАН. В качестве наблюдений рассматривались данные о полном интегральном электронном содержании вдоль заданных траекторий. Сформулирована общая постановка задачи в дифференциальной форме, исследована ее разрешимость. На основе регуляризованной постановки построен итерационный алгоритм решения задачи ассимиляции, показана его сходимость. Построена конечномерная аппроксимация и реализовано численное решение задачи, доказаны устойчивость и сходимость разностной схемы. На основе контрольных численных экспериментов исследовано качество восстановления полей распределения электронной концентрации, показано приемлемое восстановление слабовозмущенного решения как для стационарной, так и для эволюционной постановок при вертикальных и наклонных траекториях интегрирования. Библ. 38. Фиг. 8.

Ключевые слова: ионосфера, амбиполярная диффузия, обратные задачи, вариационное усвоение данных, численное моделирование.

1. ВВЕДЕНИЕ

Известно, что знание основных характеристик земной ионосферы имеет существенное практическое значение, поскольку они определяют условия для распространения всех типов радиосигналов в околоземной среде, обеспечивающих бесперебойную работу систем дальней глобальной радиосвязи, спутниковых систем связи, радиолокации, а также навигационных систем глобального спутникового позиционирования [1]–[3]. Таким образом, по мере значительного роста использования систем навигации и спутниковой радиосвязи в социальной, экономической и оборонной сферах в последние десятилетия все острее признается потребность в надежной и точной системе анализа и прогнозирования состояния ионосферы и его изменчивости. Недостаточный на сегодняшний день уровень точности описания физических процессов и высокой изменчивости в верхней атмосфере и ионосфере Земли представляет собой фундаментальную сложность решения данной задачи. Результаты последних исследований показывают, что многие механизмы формирования изменчивости остаются совсем не изученными [4], [5], а для их детального описания, по-видимому, требуется комплексный подход, учитывающий внешние воздействия космической природы (солнечная активность и др.), электромагнитные и динамические процессы во всей атмосфере, а также нелинейное взаимодействие термосферы и ионосферы [6]–[8].

Поскольку в последние десятилетия произошло существенное увеличение объема и доступности глобальных данных наблюдений о параметрах ионосферы различного типа, это сильно стимулировало различные научные группы перейти к созданию систем усвоения данных для ионосферы, а также к новому направлению по включению верхних слоев атмосферы и ионосферы в разрабатывающиеся модели Земной системы. Для решения таких задач в Институте вычислительной математики им. Г.И. Марчука (ИВМ РАН) [9] в течение последних лет уже разработаны новые версии глобальной модели общей циркуляции атмосферы (0–130 км) с включением нижних слоев термосферы и ионосферы [10–12], а также новые модели термосферы (90–500 км) [13] и F слоя ионосферы [14]. На их основе создана новая глобальная динамическая совместная модель термосферы и ионосферы [15].

Предлагаемая работа является первой из серии работ, посвященных обозначенной в начале задаче разработки перспективной прогностической системы четырехмерного вариационного усвоения данных для модели ионосферы Земли ИВМ РАН для решения проблемы прогноза глобального состояния ионосферы. Как известно, системы усвоения данных позволяют подготовить согласованные трехмерные распределения характеристик среды и начальное состояние, необходимое при решении задачи краткосрочного прогноза ее состояния. В силу небольшого времени предсказуемости первого рода для ионосферной динамики в F слое времена краткосрочного прогноза лежат в рамках суточного цикла (практическую значимость имеют времена от 15 мин до 4–6 ч). Успешное решение задачи ассимиляции данных для описания глобального состояния ионосферы стало возможным благодаря появлению и расширению ряда систем непрерывного мониторинга ее характеристик, отвечающих требованиям достаточной точности и глобальности покрытия, оперативности и доступности. В целом можно выделить два основных класса данных наблюдений ионосферы, отвечающих приведенным требованиям:

• Данные о полном содержании электронов (ПЭС/TEC) вдоль наклонных трасс связи наземных станций и спутников в фиксированные моменты времени (интеграл электронной плотности вдоль трассы радиосигнала станция-спутник), полученные с помощью обработки радиосигналов глобальных навигационных спутниковых систем (ГНСС/GNSS) или других спутниковых систем [16] (в первую очередь это данные систем GPS [16] и ГЛОНАСС [18]. Обе системы представляют собой большую группу космических аппаратов, обращающихся на среднем расстоянии от поверхности Земли 20 000 км, наземный сегмент содержит широкую сеть стационарных приемников, расположенных в основном на суше. Станции принимают спутниковый сигнал и публикуют протоколы обмена данными (с частотой обновления 30 с), по которым можно определить ПЭС вдоль траектории сигнала с точностью до некоторой неопределенной аддитивной постоянной. Для исключения этой постоянной используются отдельные методы [19]–[21], которые в данной работе обсуждаться не будут. Таким образом, в общем доступе имеется оперативная непрерывная информация о наклонных значениях ПЭС с рядом ограничений, однако с широким покрытием и высокой частотой получения, которая с использованием дополнительных методов обработки и интерполяции может служить источником информации о глобальном состоянии ионосферы.

• Данные о профилях электронной плотности, главным образом полученные ионозондами, работающими в реальном времени и количество которых в мире непрерывно увеличивается [22], [23]. C помощью ионозондов на основе отражения радиосигнала удается в реальном времени (с частотой порядка 15 минут) восстанавливать вертикальный профиль электронной плотности вплоть до высоты ее максимума в месте расположения прибора, однако их возможности ограничены низкой автоматизацией обработки данных, неточностью при низких значениях концентрации в нижних слоях ионосферы.

В разрабатываемой в ИВМ РАН системе усвоения данных предполагается использование двух описанных видов данных непрерывного мониторинга ионосферы.

Разработка систем ассимиляции данных для ионосферы была начата за рубежом в конце 1990-х годов в США в рамках программы MURI с разработки двух разных ассимиляционных численных моделей ионосферы – Global Assimilation of Ionospheric Measurements с участием групп из Университета штата Юта (USU GAIM) [24] и Университета Южной Калифорнии (USC GAIM) [25], использующих методы фильтров Калмана для усвоения данных в физические модели ионосферы. В дальнейшем стали использовать ансамблевый подход и более развитые методы. В настоящее время версия модели USU GAIM на основе ансамблевого фильтра Гаусса-Маркова используется метеорологическим агентством ВВС США (AFWA). Несколько оперативных центров космической погоды внедрили методы усвоения данных для спецификации и прогноза характеристик ионосферы на основе ассимиляции информации в эмпирические модели: ряд центров использует сервис построения оперативных карт TEC Австралийской службы прогнозов ионосферы (IPS)22, другими аналогичными примерами являются системы MIDAS-GPS [26] и EDAM-GPS [27]. Существует также ряд региональных систем усвоения данных, в основном созданных для объективного анализа и мониторинга ПЭС (к примеру US-TEC и др. [28]). В России попытка разработки ансамблевой системы усвоения данных для ионосферы на основе фильтра Калмана и американской модели CTIP была предпринята в Центральной Аэрологической Обсерватории (ЦАО) Росгидромета [29].

Особенностью подхода к созданию системы усвоения данных для модели ионосферы, излагаемой в данной работе, являются использование четырехмерной вариационной ассимиляции данных наблюдений на основе оригинальной глобальной динамической модели F слоя ионосферы [14], [30] и ее детальное математическое исследование. В качестве первой версии системы усвоения рассматривается задача ассимиляции данных на основе системы уравнений амбиполярной диффузии в квазидвумерной постановке [14]. Подобная приближенная модель достаточно хорошо описывает состояние ионосферы средних и низких широт в спокойных условиях [3], [13], [14]. Данная постановка задачи усвоения имеет смысл, поскольку в модели F слоя ионосферы ИВМ РАН для интегрирования по времени используется метод расщепления по физическим процессам [15], позволяющий включить усвоение на одном из этапов расщепления в рамках реализации всей модели. Построение системы усвоения данных на основе такого подхода для полной трехмерной модели ионосферы является предметом следующих работ. Отметим, что в простейшей версии трехмерной модели, не включающей полного описания трехмерного переноса, можно использовать предложенный в данной статье подход практически без изменений, поскольку рассматриваемые уравнения записываются одинаково на всех долготах, а зависимость от долготы и времени входит в коэффициенты уравнений, их правые части, а также в данные наблюдений. В перспективе подобный подход позволит построить системы усвоения данных для совместных моделей термосферы – ионосферы.

Основная цель данной работы – разработка и реализация первой версии системы вариационного усвоения данных для диффузионной модели ионосферы [14]. Для получающейся при этом системы уравнений требуется построить ее аппроксимацию, исследовать разрешимость и сходимость, предложить алгоритм реализации, а также провести контрольные численные эксперименты для исследования точности предложенного алгоритма. В качестве данных наблюдений в работе рассматриваются данные о ПЭС.

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

2. ПОСТАНОВКА ЗАДАЧИ

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

Пусть $z$ – высота, $\varphi $ – широта, $t$ – время, $\Omega = \{ (z,\varphi )\,{\text{|}}\,{{z}_{b}} \leqslant z \leqslant {{z}_{t}};$ $ - \pi {\text{/}}2 \leqslant \varphi \leqslant \pi {\text{/}}2\} \subset {{\mathbb{R}}^{2}}$ – сферический слой, ${{Q}_{T}} = (0,T) \times \Omega $. Рассмотрим уравнение амбиполярной диффузии в приближении тонкого сферического слоя:

$\begin{gathered} \frac{{\partial n}}{{\partial t}} = \frac{1}{{{{a}^{2}}\cos \varphi }}\frac{\partial }{{\partial \varphi }}\left[ {D{{{\cos }}^{2}}I\frac{{\partial n}}{{\partial \varphi }}\cos \varphi } \right] + \frac{\partial }{{\partial z}}\left[ {D{{{\sin }}^{2}}I\frac{{\partial n}}{{\partial z}}} \right] - \frac{1}{{a\cos \varphi }}\frac{\partial }{{\partial \varphi }}\left[ {D\sin I\cos I\frac{{\partial n}}{{\partial z}}\cos \varphi } \right] - \\ - \;\frac{1}{a}\frac{\partial }{{\partial z}}\left[ {D\sin I\cos I\frac{{\partial n}}{{\partial \varphi }}} \right] - \frac{1}{{a\cos \varphi }}\frac{\partial }{{\partial \varphi }}\left[ {\frac{D}{H}\sin I\cos I{{n}_{i}}\cos \varphi } \right] + \frac{\partial }{{\partial z}}\left[ {\frac{D}{H}{{{\sin }}^{2}}I{{n}_{i}}} \right] + P - kn, \\ \end{gathered} $
$\begin{gathered} n{{{\text{|}}}_{{t = 0}}} = {{n}_{0}}, \\ {{\left. {\left( {D{{{\sin }}^{2}}I\frac{{\partial n}}{{\partial z}} - \frac{1}{a}D\sin I\cos I\frac{{\partial n}}{{\partial \varphi }} + \frac{D}{H}{{{\sin }}^{2}}In} \right)} \right|}_{{z = {{z}_{b}}}}} = {{F}_{{lb}}}, \\ {{\left. {\left( {D{{{\sin }}^{2}}I\frac{{\partial n}}{{\partial z}} - \frac{1}{a}D\sin I\cos I\frac{{\partial n}}{{\partial \varphi }} + \frac{D}{H}{{{\sin }}^{2}}In} \right)} \right|}_{{z = {{z}_{t}}}}} = {{F}_{{ub}}}. \\ \end{gathered} $
Здесь $n(z,\varphi ,t)$ – электронная концентрация, $I = (2\varphi )$ – угол магнитного наклонения, $P(z,\varphi ,t)$ – функция фотоионизации, $k(z,\varphi ,t)$ – функция рекомбинации, ${{F}_{{ub}}}(\varphi ,t)$ и ${{F}_{{lb}}}(\varphi ,t)$ – потоки на верхней и нижней границах расчетной области, которые в данной работе мы будем полагать равными нулю. На полюсах (при $\varphi = \pm \frac{\pi }{2}$) эффективные коэффициенты диффузии равны нулю, и в задаче ставится условие ограниченности производных $\frac{{\partial n}}{{\partial \varphi }}$ по широте.

Будем считать правую часть $P$ равной ${{P}_{0}} + U$, где ${{P}_{0}}$ – известная функция фотоионизации, а $U$ – добавочное слагаемое, рассматриваемое в данной постановке задачи в качестве управления. Считаем также, что заданы данные наблюдений – конечный набор интегралов от решения вдоль заданных прямолинейных траекторий ${{\ell }_{k}}$, соединяющих спутник с точкой на поверхности Земли (TEC – total electron content или ПЭС – полное электронное содержание). Для простоты будем вести интегрирование по тонким прямоугольникам ${{\Omega }_{k}}$ (аппроксимирующим прямые ${{\ell }_{k}}$, вдоль которых рассчитывается сигнал для вычисления TEC-ов):

(1)
$\int\limits_{{{\Omega }_{k}}} n(z,\varphi ,t)a\cos \varphi d\varphi dz = Te{{c}_{k}}(t)\quad \forall t \in (0,T),\quad k = 1, \ldots ,N.$
Значения интегралов с течением времени меняются, поскольку функция $n(z,\varphi ,t)$ зависит от времени, поэтому $Te{{c}_{k}}(t)$ также зависят от времени. Эти зависимости мы полагаем известными.

Рассмотрим задачу восстановления профилей по наблюдаемым ПЭС следующего вида: найти функцию $n(z,\varphi ,t)$ и управление $U(z,\varphi ,t)$, определенные в ${{Q}_{T}}$, такие, что $n$ является решением уравнения амбиполярной диффузии c известным начальным условием и удовлетворяет данным наблюдений (1): известны интегралы по подобластям ${{\Omega }_{k}}$ и их зависимости от времени (TEC-и) (этих подобластей всего $N$).

Введем обозначения $K_{1}^{2} = D{{\sin }^{2}}I$, $K_{2}^{2} = D{{\cos }^{2}}I{{\cos }^{2}}\varphi $, $K = \left( {\begin{array}{*{20}{c}} {K_{1}^{2}}&{}&{ - {{K}_{1}}{{K}_{2}}} \\ { - {{K}_{1}}{{K}_{2}}}&{}&{K_{2}^{2}} \end{array}} \right)$ – симметричная матрица коэффициентов диффузии, а также оператор $\frac{\partial }{{\partial y}} = \frac{1}{{a\cos \varphi }}\frac{\partial }{{\partial \varphi }}$ и градиент $\nabla = \left( {\begin{array}{*{20}{c}} {\partial {\text{/}}\partial z} \\ {\partial {\text{/}}\partial y} \end{array}} \right)$. Обозначим также вертикальные и горизонтальные скорости $u = \frac{D}{H}{{\sin }^{2}}I$, ${v} = \frac{D}{H}\sin I\cos I$. Введем также для краткости обозначение $d\Omega = a\cos \varphi d\varphi dz$ – элемент площади. Тогда выписанную систему можно записать более кратко:

$\begin{gathered} \left( {\frac{\partial }{{\partial t}} - \nabla (K\nabla ({\kern 1pt} \cdot {\kern 1pt} )) - \frac{\partial }{{\partial z}}(u \cdot ({\kern 1pt} \cdot {\kern 1pt} )) - \frac{\partial }{{\partial y}}({v}\cos \varphi \cdot ({\kern 1pt} \cdot {\kern 1pt} )) + k} \right)n = {{P}_{0}} + U, \\ {{\left. n \right|}_{{t = 0}}} = {{n}_{0}},\quad {{\left. {\left( {K_{1}^{2}\frac{{\partial n}}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial n}}{{\partial y}} - un} \right)} \right|}_{{z = {{z}_{b}}}}} = 0,\quad {{\left. {\left( {K_{1}^{2}\frac{{\partial n}}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial n}}{{\partial y}} - un} \right)} \right|}_{{z = {{z}_{t}}}}} = 0, \\ \int\limits_{{{\Omega }_{k}}} n(z,\varphi ,t)d\Omega = Te{{c}_{k}}(t)\quad \forall t \in (0,T),\quad k = 1, \ldots ,N. \\ \end{gathered} $

Уравнение амбиполярной диффузии-плазмохимии для модели ионосферы подробно рассмотрено в работе [14]. Отметим ряд его ключевых особенностей.

1. Величина коэффициента диффузии $D$ меняется на 6–7 порядков при изменении высоты $z$ от $100$ до $500$ км.

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

3. Диффузия, описываемая этим уравнением, идет вдоль магнитных силовых линий. Это утверждение эквивалентно вырожденности в каждой точке тензора “коэффициентов диффузии”

$\left( {\begin{array}{*{20}{c}} {K_{1}^{2}}&{ - {{K}_{1}}{{K}_{2}}} \\ { - {{K}_{1}}{{K}_{2}}}&{K_{2}^{2}} \end{array}} \right).$
Умножив уравнение скалярно на $n$ и проинтегрировав по всей расчетной области, можно получить (при однородных краевых условиях и нулевых функциях $P$ и $k$) интегральное соотношение, отражающее геометрические характеристики диффузионных процессов (под знаком двойного интеграла в правой части стоит квадрат производной электронной концентрации по направлению магнитной силовой линии):
(2)
$\frac{1}{2}\frac{\partial }{{\partial t}}\iint\limits_{\varphi ,z} {{{n}^{2}}\cos \varphi d\varphi dz} = - \iint\limits_{\varphi ,z} {{{{\left( {{{K}_{1}}\frac{{\partial n}}{{\partial z}} - \frac{{{{K}_{2}}}}{{a\cos \varphi }}\frac{{\partial n}}{{\partial \varphi }}} \right)}}^{2}}}\cos \varphi d\varphi dz \leqslant 0.$
Корректная аппроксимация диффузионного оператора, обладающая разностным аналогом этого соотношения (и потому абсолютно устойчивая при подходящей аппроксимации по времени), уже изучена в работе [31]. Эта аппроксимация в данной работе будет использована при численном решении системы прямых и сопряженных уравнений.

3. ОБОБЩЕННАЯ ПОСТАНОВКА

Будем рассматривать функции ${{P}_{0}},U \in {{\mathbb{L}}_{2}}({{Q}_{T}})$, $Te{{c}_{k}}(t) \in {{\mathbb{L}}_{2}}(0,T)$, а решения $n$ из подпространства $W$ в пространстве ${{\mathbb{L}}_{2}}(0,T;\mathbb{W}_{2}^{1}(\Omega ))$, где $W$ состоит из таких функций, что $\frac{\partial }{{\partial t}}n \in $ $ \in {{\mathbb{L}}_{2}}\left( {0,T;(\mathbb{W}_{2}^{1}(\Omega )){\kern 1pt} *} \right)$. Производная $\frac{\partial }{{\partial t}}n$ понимается в обобщенном смысле: для всякой $n \in {{\mathbb{L}}_{2}}(0,T;\mathbb{W}_{2}^{1}(\Omega ))$ ее производная по времени есть элемент $D{\kern 1pt} '((0,T);\mathbb{W}_{2}^{1}(\Omega ))$, т.е. непрерывный оператор из пространства пробных (бесконечно гладких финитных) $\varphi \in D(0,T)$ в пространство $\mathbb{W}_{2}^{1}(\Omega )$. Его действие задается формулой $\left( {\frac{{\partial n}}{{\partial t}}} \right)(\varphi ) = - {{\left\langle {n,\left( {\frac{{\partial \varphi }}{{\partial t}}} \right)} \right\rangle }_{{{{\mathbb{L}}_{2}}[0,T]}}}$. Ограничение $\frac{\partial }{{\partial t}}n \in $ $ \in {{\mathbb{L}}_{2}}\left( {0,T;(\mathbb{W}_{2}^{1}(\Omega )){\kern 1pt} *} \right)$ означает, что требуется, чтобы производная $\frac{\partial }{{\partial t}}n$ являлась функционалом из $(\mathbb{W}_{2}^{1}(\Omega )){\kern 1pt} *$ при почти всех $t$, и нормы этих функционалов (зависящие от параметра $t$) были интегрируемы с квадратом на $[0,T]$. Известно (см. [32], гл. 3, теорема 1.1), что все функции из $W$ являются непрерывными как отображения $[0,T] \to {{\mathbb{L}}_{2}}(\Omega )$ при должном изменении на множестве меры нуль в отрезке $[0,T]$ (гарантируется, что такое изменение возможно).

Умножим скалярно первое уравнение на пробную функцию $\psi $ из класса $\mathbb{W}_{2}^{1}(\Omega )$ и проинтегрируем по $\Omega $. С помощью интегрирования по частям получим:

$\begin{gathered} \int\limits_\Omega \frac{{\partial n}}{{\partial t}}\psi d\Omega + \int\limits_\Omega \left( {K_{1}^{2}\frac{{\partial n}}{{\partial z}}\frac{{\partial \psi }}{{\partial {{z}_{{\kern 1pt} }}}}} \right. + K_{2}^{2}\frac{1}{{a\cos \varphi }}\frac{{\partial n}}{{\partial \varphi }}\frac{1}{{a\cos \varphi }}\frac{{\partial \psi }}{{\partial \varphi }} - {{K}_{1}}{{K}_{2}}\frac{1}{{a\cos \varphi }}\frac{{\partial n}}{{\partial \varphi }}\frac{{\partial \psi }}{{\partial z}} - \\ - \;\left. {{{K}_{1}}{{K}_{2}}\frac{1}{{a\cos \varphi }}\frac{{\partial n}}{{\partial z}}\frac{{\partial \psi }}{{\partial \varphi }}} \right)d\Omega + \int\limits_\Omega \,un\frac{{\partial \psi }}{{\partial z}}d\Omega + \int\limits_\Omega \,{v}n\frac{1}{a}\frac{{\partial \psi }}{{\partial \varphi }}d\Omega - \int\limits_\Omega \,(P - kn)\psi d\Omega - \\ \end{gathered} $
$\begin{gathered} \, - {{\left. {\int\limits_{ - \pi /2}^{ + \pi /2} \left( {K_{1}^{2}\frac{{\partial n}}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial n}}{{\partial y}} - un} \right)} \right|}_{{z = {{z}_{t}}}}}\psi ({{z}_{t}},\varphi ,t)a\cos \varphi d\varphi + \\ + \;{{\left. {\int\limits_{ - \pi /2}^{ + \pi /2} \left( {K_{1}^{2}\frac{{\partial n}}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial n}}{{\partial y}} - un} \right)} \right|}_{{z = {{z}_{b}}}}}\psi ({{z}_{b}},\varphi ,t)a\cos \varphi d\varphi = 0. \\ \end{gathered} $
Последнее равенство, записанное для произвольной функции $\psi \in \mathbb{W}_{2}^{1}(\Omega )$, будет обобщенной постановкой исходной задачи, если подставить $0$ вместо потоков
$\left( {K_{1}^{2}\frac{{\partial n}}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial n}}{{\partial y}} - un} \right)$
при $z = {{z}_{t}}$ и $z = {{z}_{b}}$: это исключит слагаемые с интегралами по широте на нижней и верхней границах.

Обозначим

$\begin{gathered} {\text{Diff}}(n,\psi ) = \int\limits_\Omega \left( {K_{1}^{2}\frac{{\partial n}}{{\partial z}}\frac{{\partial \psi }}{{\partial z}} + K_{2}^{2}\frac{1}{{a\cos \varphi }}\frac{{\partial n}}{{\partial \varphi }}\frac{1}{{a\cos \varphi }}\frac{{\partial \psi }}{{\partial \varphi }}} \right. - \\ \, - {{K}_{1}}{{K}_{2}}\frac{1}{{a\cos \varphi }}\frac{{\partial n}}{{\partial \varphi }}\frac{{\partial \psi }}{{\partial z}} - \left. {{{K}_{1}}{{K}_{2}}\frac{1}{{a\cos \varphi }}\frac{{\partial n}}{{\partial z}}\frac{{\partial \psi }}{{\partial \varphi }}} \right)d\Omega . \\ \end{gathered} $
Введем следующие билинейные формы и операторы:

• $\mathcal{B}(n,\psi ) = {\text{Diff}}(n,\psi ) + {{\left\langle {un,\frac{{\partial \psi }}{{\partial z}}} \right\rangle }_{{{{\mathbb{L}}_{2}}(\Omega )}}} + {{\left\langle {\frac{1}{a}{v}n,\frac{{\partial \psi }}{{\partial \varphi }}} \right\rangle }_{{{{\mathbb{L}}_{2}}(\Omega )}}} + {{\left\langle {kn,\psi } \right\rangle }_{{{{\mathbb{L}}_{2}}(\Omega )}}}$.

Определим оператор $B$ равенством $\mathcal{B}(n,\psi ) = \langle Bn,\psi \rangle $, т.е. $Bn$ – функционал, действующий на пробную функцию $\psi $. Оператор $B$ действует из пространства ${{\mathbb{L}}_{2}}\left( {0,T;(\mathbb{W}_{2}^{1}(\Omega ))} \right)$ в ${{\mathbb{L}}_{2}}\left( {0,T;(\mathbb{W}_{2}^{1}(\Omega )){\kern 1pt} *} \right)$, его область определения есть $D(B) = W$;

Наконец, определим оператор $L = \frac{\partial }{{\partial t}} + B:{{\mathbb{L}}_{2}}\left( {0,T;(\mathbb{W}_{2}^{1}(\Omega ))} \right) \to {{\mathbb{L}}_{2}}\left( {0,T;(\mathbb{W}_{2}^{1}(\Omega )){\kern 1pt} *} \right),$ $D(L) = W.$

Отметим, что повторное интегрирование по частям в обобщенной постановке приводит к соотношениям, которые впоследствии будут использованы в качестве краевых условий сопряженной задачи (они имеют такой же вид, как краевые условия для $n$ в прямой задаче, но без слагаемого $( - un)$):

$K_{1}^{2}\frac{{\partial \psi }}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial \psi }}{{\partial y}} = K_{1}^{2}\frac{{\partial \psi }}{{\partial z}} - \frac{{{{K}_{1}}{{K}_{2}}}}{{a\cos \varphi }}\frac{{\partial \psi }}{{\partial \varphi }} = 0.$
Физический смысл этого краевого условия тот же: поток вдоль магнитной силовой линии ра-вен нулю.

Сопряженный оператор

$L{\kern 1pt} * = - \frac{\partial }{{\partial t}} - \nabla (K\nabla ({\kern 1pt} \cdot {\kern 1pt} )) + u\frac{\partial }{{\partial z}} + \frac{1}{a}v\frac{\partial }{{\partial \varphi }} + k$
определен на подпространстве

$D(L{\kern 1pt} *) \subset \left\{ {q{{{\left| {{{{\left. q \right|}}_{{t = T}}} = 0,\;K_{1}^{2}\frac{{\partial \psi }}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial \psi }}{{\partial y}}} \right|}}_{{z = {{z}_{b}},z = {{z}_{t}}}}} = 0} \right\}$.

• Оператор ${{C}_{{{\text{tec}}}}}:{{\mathbb{L}}_{2}}({{Q}_{T}}) \to {{({{\mathbb{L}}_{2}}(0,T))}^{N}}$, действие которого определяется формулой ${{C}_{{{\text{tec}}}}}n = $ $ = {\mathbf{Tec}}(t)$, где ${\mathbf{Tec}}(t)$ – вектор-функция, компоненты которой есть $Te{{c}_{k}}(t)$, $k = 1,N$. Из билинейной формы

${{\left\langle {{{C}_{{{\text{tec}}}}}n,{\mathbf{f}}} \right\rangle }_{{{{{({{\mathbb{L}}_{2}}(0,T))}}^{N}}}}} = \int\limits_0^T \left( {\begin{array}{*{20}{c}} {\int\limits_{{{\Omega }_{1}}} nd\Omega } \\ \ldots \\ {\int\limits_{{{\Omega }_{N}}} nd\Omega } \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{{f}_{1}}} \\ \ldots \\ {{{f}_{N}}} \end{array}} \right)dt = \sum\limits_{k = 1}^N {\kern 1pt} {\kern 1pt} \int\limits_0^T {\kern 1pt} \int\limits_{{{\Omega }_{k}}} {\kern 1pt} n{{f}_{k}}d\Omega dt = \sum\limits_{k = 1}^N {\kern 1pt} {\kern 1pt} \int\limits_{{{Q}_{T}}} {\kern 1pt} n(z,\varphi ,t){{f}_{k}}(t){{\chi }_{{{{\Omega }_{k}}}}}(z,\varphi )d{{\Omega }_{T}}$
легко видеть, что действие оператора $C_{{{\text{tec}}}}^{*}:({{\mathbb{L}}_{2}}(0,T{{))}^{N}} \to {{\mathbb{L}}_{2}}({{Q}_{T}})$ определяется формулой:
$C_{{{\text{tec}}}}^{*}{\mathbf{f}} = \sum\limits_{k = 1}^N \,{{f}_{k}}(t){{\chi }_{{{{\Omega }_{k}}}}}(z,\varphi ),$
где ${{\chi }_{{{{\Omega }_{k}}}}}(z,\varphi )$ есть характеристическая функция подобласти ${{\Omega }_{k}}$.

Мы будем рассматривать действие ${{C}_{{{\text{tec}}}}}$ лишь на функции из пространства $W$, а образы элементов ${{({{\mathbb{L}}_{2}}(0,T))}^{N}}$ под действием $C_{{{\text{tec}}}}^{*}$ рассматривать в том числе как элементы $W{\kern 1pt} *$.

Исходную задачу теперь можно сформулировать в виде

$\begin{gathered} \langle Ln,\psi \rangle = {{\langle {{P}_{0}} + U,\psi \rangle }_{{{{\mathbb{L}}_{2}}(\Omega )}}}\quad \forall \psi \in W_{2}^{1}(\Omega )\quad {\text{для}}\;{\text{п}}{\text{.в}}{\text{.}}\;\;t \in [0,T] \\ {{C}_{{{\text{tec}}}}}n = {\mathbf{Tec}}(t). \\ \end{gathered} $

Отметим, что функция потока ${{F}_{{ub}}}$ на верхней границе расчетной области по существу неизвестна и была положена равной нулю. Одним из путей модификации предложенной модели может быть рассмотрение вектор-функции управления ${\mathbf{V}} = \left( {\begin{array}{*{20}{c}} U \\ {{{F}_{{ub}}}} \end{array}} \right) \in {{\mathbb{L}}_{2}}({{Q}_{T}}) \times {{\mathbb{L}}_{2}}([0,T] \times \Gamma )$ с тем, чтобы использовать данные наблюдений для восстановления не только функции $U$, но и функции ${{F}_{{ub}}}$.

3.1. Плотная разрешимость задачи, регуляризация и итерационный процесс

Рассмотрим задачу в операторной форме:

(3)
$Ln = {{P}_{0}} + U;\quad {{C}_{{{\text{tec}}}}}n = {\mathbf{Tec}}(t).$

Если исключить из системы (3) переменную $n$, оставив только управление $U$, то останется задача вида $\mathcal{A}U = f$, где $\mathcal{A} = C{{L}^{{ - 1}}}$, $f = {\mathbf{Tec}} - C{{L}^{{ - 1}}}{{P}_{0}}$ (здесь и далее $C = {{C}_{{{\text{tec}}}}}$ – оператор наблюдений).

Утверждение 1. Имеет место плотная разрешимость задачи $\mathcal{A}U = f$.

Доказательство. Плотная разрешимость (плотность $\operatorname{Im} \mathcal{A}$ в пространстве, куда действует $\mathcal{A}$) эквивалентна в силу теоремы Фредгольма $\operatorname{Ker} \mathcal{A}{\kern 1pt} * = 0$.

Покажем, что уравнение $\mathcal{A}{\kern 1pt} *{\kern 1pt} w = 0$ имеет только тривиальное решение. Будем при этом считать, что функции скоростей $u$ и $v$, входящие в оператор $L$, таковы, что $L$ имеет ограниченный обратный ${{L}^{{ - 1}}}$ (при $u = v = 0$ это очевидно, поскольку функция рекомбинации в правой части прямого уравнения $k \geqslant 0$, и поэтому оператор $L$ ограничен снизу). Тогда $\mathcal{A}{\kern 1pt} * = ({{L}^{{ - 1}}}){\kern 1pt} *{\kern 1pt} C{\kern 1pt} *$, и уравнение $\mathcal{A}{\kern 1pt} *{\kern 1pt} w = ({{L}^{{ - 1}}}){\kern 1pt} *{\kern 1pt} C{\kern 1pt} *{\kern 1pt} w = 0$ после обозначения $q = ({{L}^{{ - 1}}}){\kern 1pt} *{\kern 1pt} C{\kern 1pt} *{\kern 1pt} w$ принимает форму (краевые условия не указаны)

$L{\kern 1pt} *{\kern 1pt} q = C{\kern 1pt} *{\kern 1pt} w,\quad q = 0,$
или (в развернутой форме)
$\begin{gathered} \left( { - \frac{\partial }{{\partial t}} - \nabla (K\nabla ({\kern 1pt} \cdot {\kern 1pt} )) + u\frac{\partial }{{\partial z}} + \frac{1}{a}v\frac{\partial }{{\partial \varphi }} + k} \right)q = \sum\limits_{k = 1}^N \,{{\chi }_{{{{\Omega }_{k}}}}}(z,\varphi ){{w}_{k}}(t), \\ {{\left. q \right|}_{\Omega }} = 0. \\ \end{gathered} $
Из $q{{{\text{|}}}_{\Omega }} = 0$ сразу следует, что $\sum\nolimits_{k = 1}^N \,{{\chi }_{{{{\Omega }_{k}}}}}(z,\varphi ){{w}_{k}}(t) = 0$. Будем считать, что ${{\Omega }_{k}}$ попарно различны. Тогда все функции ${{w}_{k}}(t)$ тождественно равны нулю, что и требовалось.

Замечание 1. Задача $\mathcal{A}U = f$ не является однозначно разрешимой.

В самом деле, однозначная разрешимость задачи эквивалентна $\operatorname{Ker} \mathcal{A} = 0$. Несложно заметить, что пространство $\operatorname{Ker} \mathcal{A}$ отлично от нулевого, поэтому решение задачи не является единственным.

Включим рассматриваемую задачу (3) в семейство задач оптимального управления вида

(4)
$L{{n}_{\alpha }} = {{P}_{\alpha }};\quad {{J}_{\alpha }}({{P}_{\alpha }},n({{P}_{\alpha }})) = \mathop {\inf }\limits_U {{J}_{\alpha }}(U,n(U)),$
где ${{P}_{\alpha }} = {{P}_{0}} + U$, $U$ – управление, а

${{J}_{\alpha }}(U,n(U)) = \alpha \left\| U \right\|_{{{{\mathbb{L}}_{2}}({{Q}_{T}})}}^{2} + \left\| {Cn(U) - {\mathbf{Tec}}} \right\|_{{{{\mathbb{L}}_{2}}(0,T)}}^{2} = \alpha \int\limits_0^T {\kern 1pt} \int\limits_Q {{U}^{2}}d\Omega dt + \int\limits_0^T {\kern 1pt} \sum\limits_{k = 1}^N {{\left( {\int\limits_{{{\Omega }_{k}}} nd\Omega - Te{{c}_{k}}(t)} \right)}^{2}}dt.$

Тем самым рассматривается регуляризация по Тихонову $\alpha U + \mathcal{A}{\kern 1pt} *{\kern 1pt} \mathcal{A}U = \mathcal{A}{\kern 1pt} *{\kern 1pt} f$ для исходной задачи $\mathcal{A}U = f$. При $\alpha > 0$ регуляризованная задача корректно разрешима, поскольку оператор $\alpha I + \mathcal{A}{\kern 1pt} *{\kern 1pt} \mathcal{A}$ самосопряженный, а его спектр отделен от нуля.

Если же $\alpha = 0$, то такая задача оптимального управления эквивалентна задаче $\mathcal{A}{\kern 1pt} *{\kern 1pt} \mathcal{A}U = \mathcal{A}{\kern 1pt} *{\kern 1pt} f$ (симметризации Гаусса задачи $\mathcal{A}U = f$). Известно, что симметризация Гаусса для уравнения $\mathcal{A}U = f$ приводит к уравнению, равносильному уравнению $\mathcal{A}U = Vf$, где $V$ – ортопроектор на $\overline {{\text{Im}}\mathcal{A}} $. Таким образом, переход от уравнения $\mathcal{A}U = f$ и последующее рассмотрение задач (4) связаны с последовательностью корректно разрешимых (при $\alpha > 0$) задач, не эквивалентных, вообще говоря, исходной задаче (3). При этом в предположении о компактности оператора $\mathcal{A}$ известно [34], что при $\alpha \to + 0$ решение регуляризованной системы $U$ сходится к решению задачи при $\alpha = 0$, имеющему минимальную норму, т.е. к псевдорешению симметризованного уравнения $\mathcal{A}{\kern 1pt} *{\kern 1pt} \mathcal{A}U = \mathcal{A}{\kern 1pt} *{\kern 1pt} f$ (напомним, что ${\text{Ker}}\mathcal{A}$ нетривиально). Это значит, что при достаточно малых $\alpha $ решение регуляризованной задачи оптимального управления может рассматриваться как приближенное решение симметризованного уравнения. Запишем систему вариационных уравнений, выражающую условия оптимальности функционала ${{J}_{\alpha }}$:

$L{{n}_{\alpha }} = {{P}_{\alpha }};\quad L{\kern 1pt} *{\kern 1pt} {{q}_{\alpha }} = C{\kern 1pt} *{\kern 1pt} (C{{n}_{\alpha }} - {\mathbf{Tec}});\quad \alpha {{P}_{\alpha }} + {{q}_{\alpha }} = 0.$

В “классической форме записи” эти уравнения имеют следующий вид:

$\begin{gathered} \left( {\frac{\partial }{{\partial t}} - \nabla (K\nabla ( \cdot )) - \frac{\partial }{{\partial z}}(u \cdot ({\kern 1pt} \cdot {\kern 1pt} )) - \frac{\partial }{{\partial y}}(v\cos \varphi \cdot ({\kern 1pt} \cdot {\kern 1pt} )) + k} \right){{n}_{\alpha }} = {{P}_{\alpha }}, \\ {{\left. {{{n}_{\alpha }}} \right|}_{{t = 0}}} = {{n}_{0}}, \\ \end{gathered} $
$\begin{gathered} {{\left. {\left( {K_{1}^{2}\frac{{\partial {{n}_{\alpha }}}}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial {{n}_{\alpha }}}}{{\partial y}} - u{{n}_{\alpha }}} \right)} \right|}_{{z = {{z}_{b}}}}} = 0,\quad {{\left. {\left( {K_{1}^{2}\frac{{\partial {{n}_{\alpha }}}}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial {{n}_{\alpha }}}}{{\partial y}} - u{{n}_{\alpha }}} \right)} \right|}_{{z = {{z}_{t}}}}} = 0; \\ \left( { - \frac{\partial }{{\partial t}} - \nabla (K\nabla ({\kern 1pt} \cdot {\kern 1pt} )) + u\frac{\partial }{{\partial z}} + \frac{1}{a}v\frac{\partial }{{\partial \varphi }} + k} \right){{q}_{\alpha }} = \sum\limits_{k = 1}^N \,{{\chi }_{{{{\Omega }_{k}}}}}(z,\varphi )\left( {\int\limits_{{{\Omega }_{k}}} {{n}_{\alpha }}d\Omega - Te{{c}_{k}}(t)} \right), \\ {{\left. {{{q}_{\alpha }}} \right|}_{{t = T}}} = 0, \\ \end{gathered} $
$\begin{gathered} {{\left. {\left( {K_{1}^{2}\frac{{\partial {{q}_{\alpha }}}}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial {{q}_{\alpha }}}}{{\partial y}}} \right)} \right|}_{{z = {{z}_{b}}}}} = 0,\quad {{\left. {\left( {K_{1}^{2}\frac{{\partial {{q}_{\alpha }}}}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial {{q}_{\alpha }}}}{{\partial y}}} \right)} \right|}_{{z = {{z}_{t}}}}} = 0; \\ \alpha {{P}_{\alpha }} + {{q}_{\alpha }} = 0. \\ \end{gathered} $

Система прямых и сопряженных уравнений решается итерационно: рассматривается последовательность $n_{\alpha }^{m}$, $q_{\alpha }^{m}$ и $P_{\alpha }^{m}$, удовлетворяющая системе

$\begin{gathered} \left( {\frac{\partial }{{\partial t}} - \nabla (K\nabla ({\kern 1pt} \cdot {\kern 1pt} )) - \frac{\partial }{{\partial z}}(u \cdot ({\kern 1pt} \cdot {\kern 1pt} )) - \frac{\partial }{{\partial y}}(v\cos \varphi \cdot ({\kern 1pt} \cdot {\kern 1pt} )) + k} \right)n_{\alpha }^{m} = P_{\alpha }^{m}, \\ {{\left. {n_{\alpha }^{m}} \right|}_{{t = 0}}} = {{n}_{0}}, \\ \end{gathered} $
$\begin{gathered} {{\left. {\left( {K_{1}^{2}\frac{{\partial n_{\alpha }^{m}}}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial n_{\alpha }^{m}}}{{\partial y}} - un_{\alpha }^{m}} \right)} \right|}_{{z = {{z}_{b}}}}} = 0,\quad {{\left. {\left( {K_{1}^{2}\frac{{\partial n_{\alpha }^{m}}}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial n_{\alpha }^{m}}}{{\partial y}} - un_{\alpha }^{m}} \right)} \right|}_{{z = {{z}_{t}}}}} = 0; \\ \left( { - \frac{\partial }{{\partial t}} - \nabla (K\nabla ( \cdot )) + u\frac{\partial }{{\partial z}} + \frac{1}{a}v\frac{\partial }{{\partial \varphi }} + k} \right)q_{\alpha }^{m} = \sum\limits_{k = 1}^N \,{{\chi }_{{{{\Omega }_{k}}}}}(z,\varphi )\left( {\int\limits_{{{\Omega }_{k}}} n_{\alpha }^{m}d\Omega - Te{{c}_{k}}(t)} \right), \\ {{\left. {q_{\alpha }^{m}} \right|}_{{t = T}}} = 0, \\ \end{gathered} $
$\begin{gathered} {{\left. {\left( {K_{1}^{2}\frac{{\partial q_{\alpha }^{m}}}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial q_{\alpha }^{m}}}{{\partial y}}} \right)} \right|}_{{z = {{z}_{b}}}}} = 0,\quad {{\left. {\left( {K_{1}^{2}\frac{{\partial q_{\alpha }^{m}}}{{\partial z}} - {{K}_{1}}{{K}_{2}}\frac{{\partial q_{\alpha }^{m}}}{{\partial y}}} \right)} \right|}_{{z = {{z}_{t}}}}} = 0; \\ P_{\alpha }^{{m + 1}} = P_{\alpha }^{m} - {{\tau }_{m}}(\alpha P_{\alpha }^{m} + q_{\alpha }^{m}). \\ \end{gathered} $

Отметим, что итерационные процессы для решения систем оптимальности, связанных с задачей ассимиляции данных наблюдений, подробно рассматриваются в [34].

Параметр ${{\tau }_{m}}$ в текущей постановке выбран постоянным. Относительно сходимости итерационного процесса известно следующее [33]:

Утверждение 2. При положительном параметре регуляризации $\alpha > 0$ итерационный процесс сходится при ${{\tau }_{m}} \leqslant {{\tau }_{{{\text{opt}}}}} = \frac{2}{{2\alpha + {{C}_{0}} + {{C}_{1}}}}$, где отрезок $[{{C}_{0}},{{C}_{1}}]$ содержит спектр оператора $\mathcal{A}{\kern 1pt} *{\kern 1pt} \mathcal{A}$.

В частности, метод сходится при $\tau = \frac{2}{{2\alpha + {{{\left\| \mathcal{A} \right\|}}^{2}}}}$.

4. РАЗНОСТНЫЕ АППРОКСИМАЦИИ ОПЕРАТОРОВ ЗАДАЧИ

Приведем описание используемой конечно-разностной аппроксимации системы вариационных уравнений, выражающих условие оптимальности сеточного функционала $J_{\alpha }^{h}$:

(5)
${{L}^{h}}n_{\alpha }^{h} = P_{\alpha }^{h};\quad ({{L}^{h}}){\kern 1pt} *{\kern 1pt} q_{\alpha }^{h} = ({{C}^{h}}){\kern 1pt} *{\kern 1pt} ({{C}^{h}}n_{\alpha }^{h} - {\mathbf{Te}}{{{\mathbf{c}}}^{h}});\quad \alpha P_{\alpha }^{h} + q_{\alpha }^{h} = 0.$

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

$\frac{1}{2}\frac{\partial }{{\partial t}}\iint\limits_{\varphi ,z} {{{n}^{2}}\cos \varphi d\varphi dz} = - \iint\limits_{\varphi ,z} {{{{\left( {{{K}_{1}}\frac{{\partial n}}{{\partial z}} - {{K}_{2}}\frac{{\partial n}}{{\partial \varphi }}} \right)}}^{2}}}\cos \varphi d\varphi dz \leqslant 0,$
имеющего место для уравнения амбиполярной диффузии при нулевых $P$, $k$ и скорости $u$, а также нулевых краевых условиях. Подробное описание этой схемы приведено в работе [31]. Приведем эту схему для прямой задачи (для сопряженной задачи аппроксимация диффузионного оператора проводится аналогичным образом):

$\begin{gathered} \frac{{\partial {{n}_{{ij}}}}}{{\partial t}} = \left[ {{{{(K_{1}^{2})}}_{{i + 1/2,j}}}\frac{{{{n}_{{i + 1,j}}} - {{n}_{{ij}}}}}{{\Delta {{z}^{2}}}} - {{{(K_{1}^{2})}}_{{i - 1/2,j}}}\frac{{{{n}_{{ij}}} - {{n}_{{i - 1,j}}}}}{{\Delta {{z}^{2}}}}} \right] + \\ + \;\left[ {\frac{{{{{(K_{2}^{2})}}_{{i,j + 1/2}}}}}{{{{a}^{2}}\cos {{\varphi }_{{j + 1/2}}}}}\frac{{{{n}_{{ij + 1}}} - {{n}_{{ij}}}}}{{\Delta {{\varphi }^{2}}}} - \frac{{{{{(K_{2}^{2})}}_{{i,j - 1/2}}}}}{{{{a}^{2}}\cos {{\varphi }_{{j - 1/2}}}}}\frac{{{{n}_{{ij}}} - {{n}_{{ij - 1}}}}}{{\Delta {{\varphi }^{2}}}}\frac{{\cos {{\varphi }_{{j - 1}}}}}{{\cos {{\varphi }_{j}}}}} \right] + \\ \end{gathered} $
$\begin{gathered} + \;\left[ { - \frac{{{{{({{K}_{1}})}}_{{i + 1/2,j}}}{{{({{K}_{2}})}}_{{i,j + 1/2}}}}}{{a\cos {{\varphi }_{{j + 1/2}}}}}\frac{{{{n}_{{i + 1,j}}} - {{n}_{{ij}}}}}{{\Delta \varphi \Delta z}} + \frac{{{{{({{K}_{1}})}}_{{i + 1/2,j - 1}}}{{{({{K}_{2}})}}_{{i,j - 1/2}}}}}{{a\cos {{\varphi }_{{j - 1/2}}}}}\frac{{{{n}_{{i + 1,j - 1}}} - {{n}_{{ij - 1}}}}}{{\Delta \varphi \Delta z}}\frac{{\cos {{\varphi }_{{j - 1}}}}}{{\cos {{\varphi }_{j}}}}} \right] + \\ \, + \left[ { - \frac{{{{{({{K}_{1}})}}_{{i + 1/2,j}}}{{{({{K}_{2}})}}_{{i,j + 1/2}}}}}{{a\cos {{\varphi }_{{j + 1/2}}}}}\frac{{{{n}_{{ij + 1}}} - {{n}_{{ij}}}}}{{\Delta \varphi \Delta z}} + \frac{{{{{({{K}_{1}})}}_{{i - 1/2,j}}}{{{({{K}_{2}})}}_{{i - 1,j + 1/2}}}}}{{a\cos {{\varphi }_{{j + 1/2}}}}}\frac{{{{n}_{{i - 1j + 1}}} - {{n}_{{i - 1j}}}}}{{\Delta \varphi \Delta z}}} \right]. \\ \end{gathered} $

При использовании неявной схемы аппроксимации по времени (или схемы Кранка-Николсон) полученная разностная схема для решения двумерного уравнения является абсолютно устойчивой. На разностном уровне такая схема обладает аналогом приведенного интегрального соотношения:

$\begin{gathered} \sum\limits_{i = 1}^{{{N}_{z}}} {\kern 1pt} \sum\limits_{j = 1}^{{{N}_{\varphi }}} \left[ {\frac{{\partial {{n}_{{ij}}}}}{{\partial t}}} \right]{{n}_{{ij}}}a\cos {{\varphi }_{j}}\Delta \varphi \Delta z = \\ = - \sum\limits_{i = 1}^{{{N}_{z}}} {\kern 1pt} \sum\limits_{j = 1}^{{{N}_{\varphi }}} {{\left[ {{{{({{K}_{1}})}}_{{i + 1/2,j}}}\frac{{{{n}_{{i + 1,j}}} - {{n}_{{i,j}}}}}{{\Delta z}} - \frac{{{{{({{K}_{2}})}}_{{i,j + 1/2}}}}}{{a\cos {{\varphi }_{{j + 1/2}}}}}\frac{{{{n}_{{i,j + 1}}} - {{n}_{{i,j}}}}}{{\Delta \varphi }}} \right]}^{2}}a\cos {{\varphi }_{j}}\Delta \varphi \Delta z. \\ \end{gathered} $
Наличие такого соотношения, с одной стороны, корректно отражает геометрию диффузионных процессов (которые идут вдоль магнитных силовых линий), а с другой – отражает диссипативность рассматриваемой задачи.

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

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

Аппроксимация интегралов вдоль прямых (интегралов по “тонким прямоугольникам”), фигурирующих в качестве условий наблюдения (наблюдаются ПЭС на отдельных траекториях) осуществляется следующим образом. Фиксируются точки $({{z}_{t}},{{\varphi }_{t}})$ и $({{z}_{b}},{{\varphi }_{b}})$ на верхней и нижней границах расчетной области – концы отрезка интегрирования. Этот отрезок разбивается на ${{N}_{{int}}}$ частей, после чего интегрирование ведется методом прямоугольников. Узлы для вычисления интегральной суммы, вообще говоря, не попадают в узлы сетки, поэтому используется стандартная четырехточечная интерполяция (вычисляется ячейка сетки, в которую попадает соответствующий узел).

В результате для каждой фиксированной траектории, вдоль которой вычисляется ПЭС, формируется матрица коэффициентов ${{c}_{{ij}}}$, с которыми следует вычислять линейную форму $\sum\nolimits_{i = 1}^{{{N}_{z}}} {\kern 1pt} \sum\nolimits_{j = 1}^{{{N}_{\varphi }}} \,{{c}_{{ij}}}{{n}_{{ij}}}$, аппроксимирующую интеграл $\int_{{{\Omega }_{k}}} {nd\Omega } $. При этом та же сеточная функция ${{c}_{{ij}}}$ определяет и правую часть сопряженного уравнения: от фиксированной $k$-й траектории в правую часть сопряженного уравнения пойдет функция ${{c}_{{ij}}}$, умноженная на разность приближенно вычисленного ПЭС ($\sum\nolimits_{i = 1}^{{{N}_{z}}} {\kern 1pt} \sum\nolimits_{j = 1}^{{{N}_{\varphi }}} \,{{c}_{{ij}}}{{n}_{{ij}}}$) и наблюдаемого значения ПЭС.

Приведем теперь ряд утверждений об аппроксимационных свойствах регуляризованной (т.е. при $\alpha > 0$) конечно-разностной задачи (5). Прежде всего докажем, что решение конечномерной системы (5) сходится к решению дифференциальной системы оптимальности (4). Поскольку построенные разностные схемы аппроксимируют операторы дифференциальной задачи, для этого остается доказать их устойчивость. Покажем, что она имеет место при условии $\frac{{{{T}^{2}}}}{\alpha }{{\left\| {{{C}^{h}}} \right\|}^{2}} < 1$, где $T$ – длина интервала времени, на котором решается задача усвоения. В случае, если константа $\frac{{{{T}^{2}}}}{\alpha }{{\left\| {{{C}^{h}}} \right\|}^{2}}$ окажется больше единицы, достаточно разбить отрезок времени $[0,T]$ на несколько подотрезков, на каждом из которых эта константа меньше единицы, и потому на каждом подотрезке рассмотренная схема устойчива.

Отметим, что исследования устойчивости разностных схем для систем прямых и сопряженных уравнений в различных постановках приведены в работе [35].

Утверждение 3. Разностная схема для системы оптимальности (5) устойчива при условии $\frac{{{{T}^{2}}}}{\alpha }{{\left\| {{{C}^{h}}} \right\|}^{2}} < 1$.

Доказательство приведем для случая нулевой скорости $u$ (в этом случае оператор $L$ включает только производную по времени и диффузионный оператор $B$). Обозначим через ${{B}^{h}}$ конечно-разностную аппроксимацию диффузионного оператора, соответствующую рассматриваемой схеме. Будем рассматривать сферическую сеточную норму (${{\ell }_{2}}$-норму с весом $a\cos {{\varphi }_{j}}$): в этой метрике ${{B}^{h}}$ симметричен и неотрицательно определен (см. [31]).

При нулевых начальных и краевых условиях система конечно-разностных прямых и сопряженных уравнений имеет вид:

$\begin{gathered} {{n}^{{j + 1}}} - {{n}^{j}} = - \tau {{B}^{h}}{{n}^{{j + 1}}} + \tau {{U}^{j}},\quad {{n}^{0}} = {{n}_{0}}; \\ - ({{q}^{{j + 1}}} - {{q}^{j}}) = - \tau {{B}^{h}}{{q}^{j}} + \tau ({{C}^{h}}){\kern 1pt} *{\kern 1pt} ({{C}^{h}}{{n}^{j}}),\quad {{q}^{{{{N}_{t}}}}} = 0; \\ \alpha {{U}^{j}} = - {{q}^{j}}, \\ \end{gathered} $
(аппроксимация по времени – неявная). Здесь $n$, $U$ и $q$ – векторы (их размерность равна числу узлов пространственной сетки), а верхний индекс отвечает моменту времени. Отметим, что сопряженное уравнение (относительно переменной $q$) решается в обратном направлении по времени.

После исключения $U$ с помощью последнего уравнения получаем:

$\begin{gathered} (I + \tau {{B}^{h}}){{n}^{{j + 1}}} = {{n}^{j}} - \frac{\tau }{\alpha }{{q}^{j}}; \\ (I + \tau {{B}^{h}}){{q}^{j}} = {{q}^{{j + 1}}} + \tau ({{C}^{h}}){\kern 1pt} *{\kern 1pt} ({{C}^{h}}{{n}^{j}}). \\ \end{gathered} $
Введем в рассмотрение также оператор $R = (I + \tau {{B}^{h}}{{)}^{{ - 1}}}$. Поскольку в сферической метрике ${{B}^{h}}$ обладает свойствами симметричности и неотрицательной определенности, $\left\| R \right\| \leqslant 1$. Систему можно переписать в виде
$\begin{gathered} {{n}^{{j + 1}}} = R{{n}^{j}} - \frac{\tau }{\alpha }R{{q}^{j}}; \\ {{q}^{j}} = R{{q}^{{j + 1}}} + \tau R({{C}^{h}}){\kern 1pt} *{\kern 1pt} ({{C}^{h}}{{n}^{j}}){\kern 1pt} {\kern 1pt} . \\ \end{gathered} $
Применяя рекурсивно первое уравнение системы, получаем ${{n}^{{j + 1}}} = {{R}^{{j + 1}}}{{n}^{0}} - \frac{\tau }{\alpha }\sum\nolimits_{k = 0}^j \,{{R}^{{k + 1}}}{{q}^{{j - k}}}$, откуда
$\left\| {{{n}^{{j + 1}}}} \right\| \leqslant \left\| {{{n}^{0}}} \right\| + \frac{\tau }{\alpha }\sum\limits_{k = 0}^j {{\left\| R \right\|}^{{k + 1}}}\left\| {{{q}^{{j - k}}}} \right\| \leqslant \left\| {{{n}^{0}}} \right\| + \frac{{(j + 1)\tau }}{\alpha }\mathop {\max }\limits_{k = 0,{{N}_{t}}} \left\| {{{q}^{k}}} \right\| \leqslant \left\| {{{n}^{0}}} \right\| + \frac{{{{N}_{t}}\tau }}{\alpha }\mathop {\max }\limits_{k = 0,{{N}_{t}}} \left\| {{{q}^{k}}} \right\|,$
а тогда, поскольку правая часть не зависит от $j$, справедлива и оценка
$\mathop {\max }\limits_{j = 0,{{N}_{t}}} \left\| {{{n}^{j}}} \right\| \leqslant \left\| {{{n}^{0}}} \right\| + \frac{{{{N}_{t}}\tau }}{\alpha }\mathop {\max }\limits_{k = 0,{{N}_{t}}} \left\| {{{q}^{k}}} \right\|.$
Отметим, что ${{N}_{t}} \cdot \tau = T$ – длительность интервала времени, на котором решается система (5). Аналогичное рекурсивное применение второго уравнения самому к себе и последующая оценка норм приводят к похожему соотношению:
$\left\| {{{q}^{j}}} \right\| \leqslant \left\| {{{q}^{{{{N}_{t}}}}}} \right\| + \tau ({{N}_{t}} - j)\left\| {({{C}^{h}}){\kern 1pt} *{\kern 1pt} ({{C}^{h}})} \right\|\mathop {\max }\limits_{k = 0,{{N}_{t}}} \left\| {{{n}^{k}}} \right\| \leqslant \left\| {{{q}^{{{{N}_{t}}}}}} \right\| + \tau {{N}_{t}}{{\left\| {{{C}^{h}}} \right\|}^{2}}\mathop {\max }\limits_{k = 0,{{N}_{t}}} \left\| {{{n}^{k}}} \right\|,$
откуда, поскольку правая часть не зависит от $j$,
$\mathop {\max }\limits_{j = 0,{{N}_{t}}} \left\| {{{q}^{j}}} \right\| \leqslant \left\| {{{q}^{{{{N}_{t}}}}}} \right\| + \tau {{N}_{t}}{{\left\| {{{C}^{h}}} \right\|}^{2}}\mathop {\max }\limits_{k = 0,{{N}_{t}}} \left\| {{{n}^{k}}} \right\|.$
Подставим оценку для ${{\max }_{{j = 0,{{N}_{t}}}}}\left\| {{{q}^{j}}} \right\|$ в правую часть оценки для ${{\max }_{{j = 0,{{N}_{t}}}}}\left\| {{{n}^{j}}} \right\|$ и получим

$\begin{gathered} \mathop {\max }\limits_{j = 0,{{N}_{t}}} \left\| {{{n}^{j}}} \right\| \leqslant \left\| {{{n}^{0}}} \right\| + \frac{T}{\alpha }\mathop {\max }\limits_{k = 0,{{N}_{t}}} \left\| {{{q}^{k}}} \right\| \leqslant \left\| {{{n}^{0}}} \right\| + \frac{T}{\alpha }\left\| {{{q}^{{{{N}_{t}}}}}} \right\| + \frac{{{{T}^{2}}}}{\alpha }{{\left\| {{{C}^{h}}} \right\|}^{2}}\mathop {\max }\limits_{j = 0,{{N}_{t}}} \left\| {{{n}^{j}}} \right\|\; \Rightarrow \\ \Rightarrow \;\left( {1 - \frac{{{{T}^{2}}}}{\alpha }{{{\left\| {{{C}^{h}}} \right\|}}^{2}}} \right)\mathop {\max }\limits_{j = 0,{{N}_{t}}} \left\| {{{n}^{j}}} \right\| \leqslant \left\| {{{n}^{0}}} \right\| + \frac{T}{\alpha }\left\| {{{q}^{{{{N}_{t}}}}}} \right\|. \\ \end{gathered} $

Замечание 2. Приведенное доказательство устойчивости остается справедливым и в случае, когда скорость $u$ отлична от нулевой, но бездивергентна. Однако численные эксперименты показывают, что схема остается устойчивой и при ненулевой скорости $u$, не являющейся бездивергентной.

Относительно сходимости сеточных решений в численных экспериментах при устремлении параметра регуляризации $\alpha $ к нулю справедливо следующее:

Утверждение 4. Сеточные решения $(n_{\alpha }^{h},P_{\alpha }^{h})$ при $\alpha \to + 0$ сходятся, причем $P_{\alpha }^{h}$ сходится к псевдорешению ${{P}^{{h, + }}}$ конечно-разностной задачи ${{\mathcal{A}}^{h}}{{P}^{h}} = {{f}^{h}}$.

Доказательство. Пусть $\{ {{v}_{1}}, \ldots ,{{v}_{n}}\} $ и $\{ {{w}_{1}}, \ldots ,{{w}_{n}}\} $ – сингулярные ортонормированные базисы оператора ${{\mathcal{A}}^{h}}$ (ортонормированные базисы из собственных векторов операторов ${{\mathcal{A}}^{h}}({{\mathcal{A}}^{h}}){\kern 1pt} *$ и $({{\mathcal{A}}^{h}}){\kern 1pt} *{\kern 1pt} {{\mathcal{A}}^{h}}$), $\sigma _{k}^{2}$ – квадраты его сингулярных чисел, ${{\sigma }_{0}} = 0$ – нулевое сингулярное число (оно присутствует, поскольку ${{\mathcal{A}}^{h}}$ имеет ненулевое ядро). Решение регуляризованного уравнения имеет вид $P_{\alpha }^{h} = (\alpha I + ({{\mathcal{A}}^{h}}){\kern 1pt} *{\kern 1pt} {{\mathcal{A}}^{h}}{{)}^{{ - 1}}}({{\mathcal{A}}^{h}}){\kern 1pt} *{\kern 1pt} {{f}^{h}} = \sum\nolimits_{k = 1}^m \frac{{{{\sigma }_{k}}\langle {{f}^{h}},{{w}_{k}}\rangle }}{{\alpha + \sigma _{k}^{2}}}{{v}_{k}}$ (заметим, что оно не содержит компонент по векторам ${{v}_{k}}$, отвечающим ${{\sigma }_{0}} = 0$). При $\alpha \to + 0$ получаем в пределе сеточную функцию $\sum\nolimits_{k = 1}^m \frac{1}{{{{\sigma }_{k}}}}\langle {{f}^{h}},{{w}_{k}}\rangle {{v}_{k}}$, совпадающую с псевдорешением ${{P}^{{h, + }}}$ задачи ${{\mathcal{A}}^{h}}{{P}^{h}} = {{f}^{h}}$, т.е. с решением задачи $({{\mathcal{A}}^{h}}){\kern 1pt} *{\kern 1pt} {{\mathcal{A}}^{h}}{{P}^{h}} = ({{\mathcal{A}}^{h}}){\kern 1pt} *{\kern 1pt} {{f}^{h}}$ наименьшей нормы.

На практике сходимость итерационного процесса для конечно-разностной задачи наблюдается и при $\alpha = 0$, но параметр ${{\tau }_{m}}$ приходится выбирать на порядок меньше, чем в регуляризованном случае, даже при достаточно малых $\alpha $.

Замечание 3. Псевдорешение (сеточная функция) ${{P}^{{h, + }}}$ аппроксимирует некоторое решение дифференциальной задачи $\mathcal{A}P = f$, при котором значение функционала ${{J}_{{\alpha = 0}}}(P)$ близко к минимальному. При этом конечно-разностная задача ${{\mathcal{A}}^{h}}{{P}^{h}} = {{f}^{h}}$ плотно разрешимой не является (ядра ${{\mathcal{A}}^{h}}$ и $({{\mathcal{A}}^{h}}){\kern 1pt} *$ ненулевые одновременно), поэтому минимальное значение сеточного функционала $J_{{\alpha = 0}}^{h}({{P}^{h}})$ отлично от нуля и равно квадрату расстояния от сеточной функции ${{f}^{h}}$ до образа $\mathcal{A}$. Заметим, однако, что при измельчении шага сетки нижняя грань множества значений $J_{\alpha }^{h}$ стремится к нулю, поскольку значение $J_{\alpha }^{h}$ при соответствующем выборе сеточных норм аппроксимирует значение ${{J}_{\alpha }}$.

5. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

Для исследования качества восстановления распределения электронной концентрации проведем серию контрольных численных экспериментов, использующих в качестве данных наблюдений тестовые значения ПЭС с заданным относительным отклонением от решения прямой задачи. Проводится модельный эксперимент по усвоению данных наблюдений следующего вида: ищется численное решение невозмущенной прямой задачи на всем отрезке времени, затем оно изменяется заданным относительно малым возмущением, после чего вычисляются интегралы (ПЭС) по выбранным траекториям. В качестве прямой модели в данной работе используем уже рассмотренную в предыдущих работах [14], [30] диффузионную модель F слоя ионосферы, все параметры которой соответствуют стандартному состоянию невозмущенной ионосферы. Тем самым мы получим тестовые данные наблюдений (на всем отрезке времени) по возмущенному распределению электронной концентрации, вообще говоря, не являющемуся решением уравнения амбиполярной диффузии. В ходе тестовых численных экспериментов будем восстанавливать с помощью системы ассимиляции решение по построенным возмущенным данным наблюдений. При этом в итерационном процессе для решения системы (5) на одной итерации решается прямое уравнение (“вперед по времени”), а затем сопряженное уравнение (“назад по времени”) на всем временном отрезке, после чего правая часть прямой задачи корректируется для перехода к следующей итерации.

Для всех численных экспериментов в качестве тестовых данных наблюдений будем использовать решения прямой задачи, возмущенные умножением на функцию $(1 + 0.1\sin (\pi (z - 100){\text{/}}400))$ (эта функция на верхней и нижней границах обращается в $1$, а внутри расчетной области меняет исходное решение не более, чем на $10\% $). В качестве траекторий интегрирования рассмотрим по отдельности случай отдельных вертикальных ПЭС и случай непрерывного заполнения части расчетной области наклонными ПЭС.

Отметим, что во всех проведенных экспериментах $10$ итераций достаточно для установления рассчитываемых профилей: на практике дальнейшие итерации не приводили к значимому изменению решения. При этом сходимость наблюдалась и при нулевом параметре регуляризации, но взятие достаточно малого положительного значения $\alpha $ позволяет на порядок увеличить параметр $\tau $ итерационного процесса без заметного изменения восстановленного численного решения.

5.1. Восстановление профилей в дневном стационарном распределении

Для исследования качества восстановления стандартного стационарного пространственного распределения электронной концентрации рассмотрим постановку прямой задачи при фиксированной правой части, для которой имеет место сходимость к стационарному решению. В этой постановке на отрезке времени 24 ч рассматривается сходимость к дневному стационарному распределению электронной концентрации с нулевых начальных условий при постоянной функции фотоионизации ${{P}_{0}}$. При такой постановке все параметры модели фиксированы во времени и соответствуют стандартному полуденному невозмущенному состоянию ионосферы в условиях весеннего равноденствия [14], [30]. Для прямой задачи при фиксированной правой части, отвечающей дневному невозмущенному состоянию системы, решение сходится к стационарному (характерное время сходимости на всей расчетной области с нулевых начальных условий – порядка $6$ ч). Применим построенный алгоритм восстановления решения по вычисленным данным наблюдения и получим решение конечно-разностной задачи, дающее приближенное решение системы оптимальности (5). Цель такого контрольного эксперимента – проверка эффективной работы алгоритма и качества восстановления пространственного распределения электронной концентрации.

Для первого эксперимента для вычислений ПЭС используем строго вертикальные траектории. Будем в качестве данных наблюдений рассматривать интегралы возмущенного решения по вертикальным траекториям на ограниченном количестве точек (на фиксированных широтах: ‒71°, –53°, –35°, –17°, 0°, 17°, 35°, 53° и 71°).

На фиг. 1 представлены результаты восстановления тестового вертикального профиля для области экватора, низких и средних широт (отметим, что при постановке предполагалось существование наблюдений на этих широтах). Эксперименты показали аналогичные результаты и для симметричных относительно экватора широт. Для всех широт интегральная ошибка много меньше (на $3$ порядка), чем интеграл от восстановленного (или от “истинного”) профиля. Отметим, что ассимиляция смещает положение максимума электронной концентрации вверх, а отклонение восстановленного профиля от “истинного” по высоте составляет $4\% $ вблизи максимума (такая ошибка по высоте расположена приблизительно на $40$ км ниже максимума) и максимально вблизи верхней границы. При этом форма профиля в целом немного изменяется, что не противоречит постановке, поскольку ядро оператора $\mathcal{A}$ ненулевое. Наибольшая ошибка восстановления вблизи верхней границы составляет около 10% и связана с выбором вида функции возмущения, на которую было умножено решение при подсчете модельных данных наблюдений: после умножения изменилась структура профилей, и возмущенное распределение электронной концентрации уже не является решением прямой задачи, а наклон профиля вблизи верхней границы не согласован с краевыми условиями для уравнения амбиполярной диффузии.

Фиг. 1.

Вертикальные распределения электронной концентрации для широт 0° (а), 35° (в) и 53° (д), полученные решением только прямой задачи для модели ионосферы (голубая линия), решением задачи восстановления (фиолетовая линия) и ожидаемое по восстановлению тестовое “точное” решение (зеленая), а также профили абсолютной разницы восстановленного и точного решения для широт 0° (б), 35° (г) и 53° (е). В качестве тестовых данных наблюдений использованы данные о вертикальных ПЭС. Концентрации электронов приведены в см–3.

На фиг. 2 приведено распределение по широтам одного из ключевых практически значимых параметров ионосферы – величины максимума (по высоте $z$) электронной концентрации и ошибки восстановления. Видно, что с точки зрения максимальной по высоте электронной концентрации предложенный алгоритм с высокой точностью (ошибка в средних широтах, на которых были выбраны данные наблюдений, не превосходит $1\% $ после ассимиляции по сравнению с исходным возмущением порядка $10\% $, а наименьшая точность $ \sim {\kern 1pt} 2\% $ получена на экваторе) восстанавливает искомое значение, причем усвоение данных по вертикали влияет на близлежащие широты приблизительно на 2° по обе стороны от широты, где расположена траектория интегрирования.

Фиг. 2.

Широтные распределения максимума электронной концентрации по высоте (а), полученные решением только прямой задачи для модели ионосферы (голубая линия), решением задачи восстановления (фиолетовая линия) и ожидаемое по восстановлению тестовое “точное” решение (зеленая), а также разница восстановленного и точного широтных распределений максимума (б). Используемые данные наблюдений – вертикальные ПЭС. Концентрации электронов приведены в см–3.

Теперь рассмотрим эксперимент с наблюдаемыми данными – интегралами вдоль заданных наклонных прямолинейных траекторий, оставив прежними все остальные параметры модели. Зададим фиксированный угол наклона прямолинейных траекторий следующим образом: точку на высоте ${{z}_{t}} = 500$ км и широте $\varphi $ соединим с точкой на высоте ${{z}_{b}} = 100$ км и широте $(\varphi - 35^\circ )$, а затем покроем максимально возможную часть области такими траекториями. При таком подходе вертикальные профили на широтах от –55° до 55° оказываются полностью покрыты траекториями интегрирования, поэтому вновь можно провести сравнение восстановленных и возмущенных профилей на этих широтах. На фиг. 3 приведены результаты восстановления профилей для тех же широт, что и в случае вертикальных профилей наблюдения.

Фиг. 3.

Вертикальные распределения электронной концентрации для широт 0° (а), 35° (в) и 53° (д), полученные решением только прямой задачи для модели ионосферы (голубая линия), решением задачи восстановления (фиолетовая линия) и ожидаемое по восстановлению тестовое “точное” решение (зеленая), а также профили абсолютной разницы восстановленного и точного решения для широт 0° (б), 35° (г) и 53° (е). В качестве тестовых данных наблюдений использованы данные о наклонных ПЭС. Концентрации электронов приведены в см–3.

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

Фиг. 4.

Широтные распределения максимума электронной концентрации по высоте (а), полученные решением только прямой задачи для модели ионосферы (голубая линия), решением задачи восстановления (фиолетовая линия) и ожидаемое по восстановлению тестовое “точное” решение (зеленая), а также разница восстановленного и точного широтных распределений максимума (б). Используемые данные наблюдений – наклонные ПЭС. Концентрации электронов приведены в см–3.

На широтах $[ - 90^\circ , - 65^\circ ]$ ассимиляция практически не повлияла на величину максимума, поскольку траектории интегрирования не затрагивали часть вертикального профиля, содержащую максимум. То же относится и к области в окрестности северного полюса. В остальных широтах восстановление максимума происходит приблизительно с одинаковой точностью (не более $3\% $ по сравнению с исходным возмущением в $10\% $). Отметим, что абсолютная ошибка восстановления максимума на экваторе втрое больше, чем в средних широтах, но и значение электронной концентрации в экваториальной зоне также приблизительно втрое больше.

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

5.2. Восстановление профилей в суточном ходе

Для более реалистичного эксперимента рассмотрим суточный ход на отрезке времени 48 ч, в котором начальные условия выбираются в соответствии с установившимся суточным циклом, фотоионизация ${{P}_{0}}$ меняется в зависимости от времени в соответствии с изменением зенитного угла Солнца [36], а остальные внешние параметры модели остаются фиксированными. Такая постановка позволяет изучить характеристики восстановления в динамике. Как и в предыдущем эксперименте, данные считаются известными на каждом шаге по времени, а одной итерации процесса решения системы прямых и сопряженных уравнений соответствует решение прямого уравнения на всем временном интервале, а затем сопряженного уравнения на том же интервале времени в обратном направлении, после чего обновляется правая часть.

Исследуем качество ассимиляции данных наблюдений при моделировании суточного хода за время $T$, равное $2$ сут. Цель эксперимента – проверка качества восстановления решения не только по пространству, но и во времени. Приведем результаты восстановления решения при задании вертикальных ПЭС на ограниченном количестве широт, аналогично предыдущим (см.  п. 5.1) экспериментам (–71°, –53°, –35°, –17°, 0°, 17°, 35°, 53° и 71°). Восстановленное решение на широтах 0°, 35° и 53°, а также разница между возмущенным и восстановленным решениями на этих широтах представлены на графиках.

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

Фиг. 5.

Временной ход вертикальных распределений электронной концентрации для широт 0° (а), 35° (в) и 53° (д), полученный решением задачи восстановления, а также профиль абсолютной разницы восстановленного и точного решения для широт 0° (б), 35° (г) и 53° (е). В качестве тестовых данных наблюдений использованы данные о вертикальных ПЭС. Полудню соответствуют отметки 12 и 36 на оси времени. Изолинии проведены в 106 и 104 см–3 соответственно.

Фиг. 6.

Временной ход вертикальных распределений электронной концентрации для широт 0° (а), 35° (в) и 53° (д), полученный решением задачи восстановления, а также профиль абсолютной разницы восстановленного и точного решения для широт 0° (б), 35° (г) и 53° (е). В качестве тестовых данных наблюдений использованы данные о наклонных ПЭС. Полудню соответствуют отметки 12 и 36 на оси времени. Изолинии проведены в 106 и 104 см–3 соответственно.

Представленные результаты показывают, что при ассимиляции данных для обоих типов ПЭС (вертикальных и наклонных) решение приближается к истинному возмущенному. Интегралы от полученного решения совпадают с выбранными данными наблюдений вплоть до третьего порядка. При этом максимум электронной концентрации немного сдвигается вверх по высоте, а наибольшая ошибка восстановления в окрестности максимума расположена приблизительно на $40$ км ниже него и не превосходит $2\% $ (на экваторе эта ошибка максимальна) по сравнению с исходной ошибкой $10\% $ для невозмущенного решения. Во всей расчетной области ошибка восстановления максимальна вблизи верхней границы, как и в стационарном случае, в силу несогласованности возмущения с верхним краевым условием прямой задачи (восстановленное решение обязано ему удовлетворять, а возмущенное – заведомо не удовлетворяет). При этом вблизи верхней границы восстановленное решение оказывается несколько завышенным. Тем самым в динамике, при переменной по времени правой части, с точки зрения характеристик восстановленных профилей наблюдаются те же эффекты, что и для восстановленного дневного стационарного решения в п. 5.1. В суточным цикле наибольшие относительные ошибки восстановления ожидаемо соответствуют дневным условиям, ошибки в ночное время не превышают 1% для всех экспериментов.

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

Фиг. 7.

Временной ход максимума распределения электронной концентрации для широт 0° (а), 35° (в) и 53° (д), полученный решением задачи восстановления для случаев вертикальных траекторий (фиолетовая линия), наклонных траекторий (зеленая линия), ожидаемая “точная” зависимость максимума от времени (голубая линия) и зависимость максимума для невозмущенного решения прямой задачи (желтая линия), а также профиль абсолютной разницы восстановленного и точного значения максимума для вертикальных траекторий (фиолетовая линия) и наклонных траекторий (зеленая линия). Представлены широты 0° (б), 35° (г) и 53° (е).

Из представленных графиков видно, что ассимиляция существенно приближает максимальное значение электронной концентрации к “истинному” (ошибка составляет не более $2\% $ для вертикальных ПЭС и $3\% $ для наклонных ПЭС соответственно по сравнению с ошибкой в $10\% $ для невозмущенного решения). При этом восстановление по наклонным траекториям оказывается несколько хуже, поскольку большая часть такой траектории расположена в окрестностях верхней и нижней границ расчетной области, в то время как максимум электронной концентрации расположен ближе к ее середине по оси $z$ (на высоте порядка 300 км).

Отметим, что более качественное восстановление решения при вертикальных ПЭС было связано в том числе с конкретным видом возмущения. В самом деле, предположим, что решение было возмущено таким образом, чтобы интеграл по любой вертикальной траектории от данного возмущения равняется нулю. Тогда в процессе усвоения по вертикальным ПЭС в правой части сопряженного уравнения стоит тождественный ноль, и начальное приближение в итерационном процессе для решения системы оптимальности (4) изменяться не будет. Тем не менее при выборе наклонных ПЭС в том же итерационном процессе решение уже будет меняться. Приведем результат восстановления максимума электронной концентрации для эксперимента, в котором данные наблюдений вычислены по решению прямой задачи $n(z,\varphi ,t)$, возмущенному по правилу ${{n}_{{{\text{new}}}}}(z,\varphi ,t) = n(z,\varphi ,t) + 0.2\sin (2\pi (z - 100){\text{/}}400)\mathop {\max }\nolimits_{z \in [100,500]} n(z,\varphi ,t)$. Нетрудно заметить, что интеграл по $z$ на отрезке от $100$ км до $500$ км от добавочного слагаемого в правой части этой формулы равен нулю. Такое возмущение изменяет форму профиля распределения электронной концентрации по высоте, причем знаки возмущения различны в областях выше и ниже высоты $z = 300$ км. Ассимиляция данных приближала решение к возмущенному в части области выше максимума электронной концентрации (ошибка с $30\% $ снижалась в этой подобласти до $15\% $ в приэкваториальной зоне, с 15 до $3\% $ вблизи широты 35° и с 15 до $7\% $ в средних широтах) и мало изменяла профили ниже $z = 300$ км, поскольку коэффициент диффузии $D$ принимает выше $z = 300$ км значения, на несколько порядков превосходящие соответствующие значения в нижней части расчетной области. Рассмотрим восстановление максимальных по $z$ значений решения на различных широтах.

На фиг. 8 зависимость максимума электронной концентрации после усвоения при выборе вертикальных траекторий интегрирования не представлена, поскольку она совпадает с невозмущенным временным ходом этой величины. При этом использование данных наблюдений на наклонных траекториях позволяет улучшить восстановление неизвестного возмущенного распределения электронной концентрации. Отметим, что выбранное возмущение уменьшает значения функции $n(z,\varphi ,t)$ в верхней части расчетной области и увеличивает в нижней, поэтому в “точном” возмущенном решении положение максимума электронной концентрации оказывается смещено по высоте. Это приводит к тому, что, в частности, максимальное значение на экваторе в возмущенном решении уменьшается, а в средних широтах увеличивается, и ошибки восстановления этой величины имеют характер, отличный от приведенного на фиг. 7 для эксперимента с ненулевым интегралом по вертикали от возмущения.

Фиг. 8.

Временной ход максимума распределения электронной концентрации для широт 0° (а), 35° (в) и 53° (д), полученный решением задачи восстановления для случая наклонных траекторий (фиолетовая линия), ожидаемая “точная” зависимость максимума от времени (зеленая линия) и зависимость максимума для невозмущенного решения прямой задачи (голубая линия), а также профиль абсолютной разницы восстановленного и точного значения максимума для наклонных траекторий (фиолетовая линия). Случай равенства нулю интеграла от возмущения решения по вертикали. Представлены широты 0° (б), 35° (г) и 53° (е).

6. ЗАКЛЮЧЕНИЕ

В заключение кратко перечислим основные результаты данной работы.

1. Сформулирована постановка задачи вариационной ассимиляции данных наблюдений для разрабатываемой в ИВМ РАН динамической модели F слоя ионосферы Земли в двумерной диффузионной постановке [14], [30]. В качестве основного типа данных в данной работе рассматриваются интегралы электронной плотности (ПЭС) вдоль прямолинейных траекторий, соединяющих спутник с точкой на земной поверхности. В качестве управления в задаче оптимальности рассматривается неопределенная добавка к правой части (функции ионизации). Приведена полная постановка задачи, показана ее плотная разрешимость и рассмотрена ее регуляризация, на основе которой построен алгоритм решения с помощью итерационного метода, доказана сходимость сеточных решений к псевдорешению конечно-разностной задачи при стремлении параметра регуляризации к нулю.

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

3. На основе контрольных численных экспериментов по восстановлению профилей по известным ПЭС для возмущений, максимум которых скоррелирован с максимумом невозмущенного распределения электронной плотности, показано высокое качество восстановления распределения электронной концентрации в областях с существованием данных наблюдений. Рассмотрены восстановление профилей в задаче о сходимости к дневному стационарному распределению электронной концентрации с нулевых начальных условий, а также восстановление профилей в суточном ходе. В обоих случаях алгоритм позволяет достаточно точно восстановить возмущенное тестовое решение (с отклонениями до $10\% $) с наибольшими ошибками 2–4% (в дневных условиях, главным образом в утренние часы) в области максимума и верхней границы (что связано с особенностями функции возмущения).

Отметим, что в данной работе рассмотрена первая версия системы усвоения данных для модели ионосферы ИВМ РАН в двумерной диффузионной постановке [14]. В качестве промежуточного варианта трехмерной модели, не использующей трехмерного переноса, может быть рассмотрен предложенный в данной работе алгоритм, в котором в качестве данных наблюдений используются интегралы по трехмерным траекториям, не локализованным на заданной долготе. На следующих этапах работы на основе разработанных алгоритмов будет реализована система усвоения для трехмерной модели ионосферы, учитывающей трехмерный перенос. Один из вариантов обобщения предложенной в статье двумерной системы усвоения на трехмерный вариант заключается в использовании разработанного алгоритма на шаге расщепления трехмерной модели. В рамках трехмерной версии системы усвоения будет отдельно рассмотрена проблема ее эффективности и возможности использования в оперативном режиме для прогноза состояния ионосферы. Отдельно в данной модели будет реализовано усвоение другого типа данных – полных или частичных вертикальных профилей (получаемых с помощью ионозондов), что потребует изменение оператора наблюдений в постановке задачи и отдельного исследования точности и эффективности такой системы при различных конфигурациях данных наблюдений.

Список литературы

  1. Дэвис К. Радиоволны в ионосфере. М.: Мир, 1973. 502 с.

  2. Иванов-Холодный Г.С., Никольский Г.М. Солнце и ионосфера: коротковолновое излучение Солнца и его воздействие на ионосферу. М.: Наука, 1969. 455 с.

  3. Kelley M.C. The Earth’s Ionosphere. Academic Press, 1989. 480 p.

  4. Deminov M.G., Deminova G.F., Zherebtsov G.A., Polekh N.M. Statistical properties of variability of the quiet ionosphere F2-layer maximum parameters over Irkutsk under low solar activity // Advances in Space Research. 2011. V. 51. № 5. P. 702–711.

  5. Ratovsky K.G., Medvedev A.V., Tolstikov M.V. Diurnal, seasonal and solar activity pattern of ionospheric variability from Irkutsk Digisonde data // Advances in Space Research. 2015. V. 55. № 8. P. 2041–2047.

  6. Miller A., Schmidt H., Bunzel F. Vertical coupling of the middle atmosphere during sudden stratospheric warming events // J. Atmos. Sol.-Terr. Phys. 2013. V. 97. P. 11–21.

  7. Yigit E., Koucki Knizova P., Georgieva K., Ward W. A review of vertical coupling in the Atmosphere–Ionosphere system: Effects of waves, sudden stratospheric warmings, space weather, and of solar activity // J. Atmos. Sol.-Terr. Phys. 2016. V. 141. P. 1–12.

  8. Pedatella N.M. et al. Multimodel comparison of the ionosphere variability during the 2009 sudden stratosphere warming // J. of Geophysical Research: Space Phys. 2016. V. 121. № 7. P. 7204–7225.

  9. Dymnikov V.P., Lykosov V.N., Volodin E. M. Mathematical simulation of earth system dynamics // Izvestiya. Atmospheric and Oceanic Phys. 2015. V. 51. № 3. P. 227–240.

  10. Kulyamin D.V., Volodin E.M. INM RAS coupled atmosphere–ionosphere general circulation model INMAIM (0–130 km) // Rus. J. of Num. Analys. and Math. Model. 2018. V. 33. № 6. P. 351–357.

  11. Кулямин Д.В., Дымников В.П. Моделирование климата нижней ионосферы // Известия Российской академии наук. Физика атмосферы и океана. 2015. Т. 51. № 3. С. 317–337.

  12. Kulyamin D.V., Dymnikov V.P. Numerical modelling of coupled neutral atmospheric general circulation and ionosphere D region // Rus. J. of Num. Analys. and Math. Model. 2016. V. 31. № 3. P. 159–171.

  13. Кулямин Д.В., Галин В.Я., Погорельцев А.И. Моделирование общей циркуляции термосферы с включением параметризации радиационных процессов // Метеорология и гидрология. 2015. Т. 40. № 6. С. 48–57.

  14. Кулямин Д.В., Останин П.А., Дымников В.П. Моделирование F слоя земной ионосферы. Решение уравнений амбиполярной диффузии // Матем. моделирование. 2019. Т. 31. № 4. С. 57–74.

  15. Дымников В.П., Кулямин Д.В., Останин П.А. Совместная модель глобальной динамики термосферы и ионосферы Земли // Известия Российской академии наук. Физика атмосферы и океана. 2020. Т. 56. № 3. С. 280–292.

  16. Hofmann-Wellenhof B., Lichtenegger H., Wasle E. Global Navigation Satellite Systems: GPS, GLONASS, Galileo and more. Wien, New York: Springer, 2008. 545 p.

  17. Parkinson B.W., Spilker J.J. Global Positioning System: Theory and Applications // Ciencia militar y naval UR. American Institute of Aeronautics & Astronautics 1996.

  18. Глобальная навигационная спутниковая система ГЛОНАСС. Интерфейсный контрольный документ. М.: 1995.

  19. Afraimovich E.L., Astafyeva E.I., Demyanov V.V., Edemskiy I.K., Gavrilyuk N.S., Ishin A.B., et al. A review of GPS/GLONASS studies of the ionospheric response to natural and anthropogenic processes and phenomena // J Space Weather Space Clim. 2013. V. 3. № A27.

  20. Rideout W., Coster A. Automated GPS processing for global total electron content data // GPS Solutions. 2006. V. 10 (3) P. 219–228. https://doi.org/10.1007/s10291-006-0029-5

  21. Brunini C., Meza A., Azpilicueta F., Zele M.A.V., Gende M., Daz A. A New Ionosphere Monitoring Technology Based on GPS // Astrophys. Space Sci. 2004. V. 290. P. 415–429.

  22. Bibl K. Evolution of the ionosonde // Annali di Geophysica. 1998. V. 41. № 5-6. P. 667–680.

  23. Reinisch B.W., Galkin I.A. Global Ionospheric Radio Observatory (GIRO) // Earth Planet Sp. 2011. V. 63. P. 377–381. https://doi.org/10.5047/eps.2011.03.001

  24. Schunk R.W., Scherliess L., Sojka J.J., Thompson D.C., Anderson D.N., Codrescu M., Minter C., Fuller-Rowell T.J., Heelis R.A., Hairston M., Howe B. Global assimilation of ionospheric measurements (GAIM) // Radio Sci. 2004. V. 39. № RS1S02.

  25. Wang C., Hajj G., Pi X., Rosen I.G., Wilson B. Development of the Global Assimilative IonosphericModel // R-adio Sci. 2004. V. 39. № RS1S06. https://doi.org/10.1029/2002RS00854

  26. Dear R., Mitchell C. GPS interfrequency biases and total electron content errors in ionospheric imaging over Europe // Radio Sci. 2006. V. 41. № RS6007. https://doi.org/10.1029/2005RS003269

  27. Angling M., Shaw J., Shukla A., Cannon P. Development of an HF selection tool based on the Electron Density Assimilative Model near-real-time ionosphere // Radio Sci. 2009. V. 44. № RS0A13. https://doi.org/10.1029/2008RS004022

  28. Fuller-Rowell T. et al. US-TEC: a new data assimilation product from the space environment center characterizing the ionospheric total electron content using real-time GPS data // Radio Sci. 2006. V. 41. https://doi.org/10.10292005RS003393

  29. Khattatov B. et al. Ionospheric nowcasting via assimilation of GPS measurements of ionospheric electron content in a global physics-based time-dependent model // Q. J. R. Meteorol. Soc. 2005. V. 131. P. 3543–3559.

  30. Ostanin P.A., Kulyamin D.V., Dymnikov V.P. Numerical modelling of the Earth’s ionosphere F region // IOP Conf. Series: Earth and Environmental Science. 2017. V. 96. № 1. P. 012011(1)–012011(11).

  31. Ostanin P.A. On the approximation of the diffusion operator in the ionosphere model with conserving the direction of geomagnetic field // Rus. J. of Num. Analys. and Math. Model. 2022. V. 37. № 1. P. 25–39.

  32. Лионс Ж.-Л. Оптимальное управление системами, описываемыми уравнениями с частными производными. М.: МИР, 1972. 416 с.

  33. Шутяев В.П. Операторы управления и итерационные алгоритмы в задачах вариационного усвоения данных. М.: Наука, 2001. 239 с.

  34. Агошков В.И. Методы оптимального управления и сопряженных уравнений в задачах математической физики. М. : ИВМ РАН, 2016. 244 с.

  35. Пармузин Е.И., Шутяев В.П. О численных алгоритмах решения одной задачи об усвоении данных // Ж. вычисл. матем. и матем. физ. 1997. V. 37. № 7. С. 816–827.

  36. Schunk R.W., Nagy A.F. IONOSPHERES Physics, Plasma Physics, and Chemistry. New York, United States: Cambridge University Press, 2009. 628 p.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал