CC BY
22
Вычислительные технологии
Том 19, № 4, 2014
Прямой метод вариационного усвоения данных
для моделей конвекции-диффузии
>к
на основе схемы расщепления*
А. В. ПЕНЕНКО, В. В. ПЕНЕНКО Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: aleks@ommgp.sscc.ru, penenko@sscc.ru
Пененко А.В., Пененко В.В. Прямой метод вариационного усвоения данных для моделей конвекции-диффузии на основе схемы расщепления // Вычислительные технологии. 2014. Т. 19, № 4. С. 69-83.
Рассматривается задача усвоения данных измерений концентраций примесей в модели конвекции-диффузии. Включение в систему моделирования данных осуществляется в рамках вариационного подхода в формулировке со слабыми ограничениями и применением схем расщепления. В модель процессов конвекции-диффузии добавляются специальные функции управления, которые находятся из условий минимума целевого функционала, объединяющего невязки между измеренными и смоделированными значениями, а также норму функций управления. Усвоение данных производится внутри стадий расщепления, что обеспечивает вычислительную эффективность схемы усвоения данных, которая не требует для своей реализации итераций при условии заданного параметра усвоения, регулирующего точность воспроизведения данных измерений.
Ключевые слова: усвоение данных, вариационный подход, конвекция-диффузия, метод расщепления.
Penenko A.V., Penenko V.V. Direct data assimilation method for convection-diffusion models based on splitting scheme // Computational Technologies. 2014. Vol. 19, No. 4. P. 69-83.
Modern urban and industry areas emit large number of airborne chemical substances that can seriously pollute the environment. These processes can be studied with multidimensional convection-diffusion models. Both high spatial dimensionality and large variety of considered chemical species impose strong requirements to the performance of direct and inverse modeling algorithms.
In this work, we consider a data assimilation problem for a nonstationary multidimensional convection-diffusion model with in situ concentration measurements. Data are assimilated at each time step of time discretized model. From theoretical point of view, data assimilation problem is considered as a sequence of inverse problems for different compositions of measurement data. Hence Tikhonov regularization can be a theoretical background to study the problem as an ill-posed measurement operator inversion problem with convection-diffusion model as a regularizer.
* Работа поддержана грантами РФФИ 14-01-31482_мол_а и 14-01-00125, Программами Президиума РАН № 4 и Отделения математических наук РАН № 3, Интеграционными проектами СО РАН № 8 и 35.
The model is augmented with control functions that are added to the source term in order to resolve uncertainty of the emission rate inherent to the atmospheric chemistry models. They are sought as the minimum of the target functional in which a control function norm is combined with the misfit between measured and corresponding modeled values.
Original element in the model lies in variational approach that has been applied to the weak-constraint formulation for each splitting stage of the multidimensional model. Because of the efficient algorithm for one-dimensional data assimilation to the convection-diffusion model, the resulting algorithm doesn't take iterations with given data assimilation parameter that regulates the measurements reproduction quality. This parameter can be chosen with Morozov discrepancy principle. Algorithm efficiency has been shown in computational experiments.
Keywords: data assimilation, variational approach, convection-diffusion, splitting method.
Введение
Цель применения алгоритмов усвоения данных — улучшение прогноза состояния некоторой системы на основе её математической модели с использованием данных измерений. Разрабатываемую технологию предполагается применять в совместных задачах динамики и химии атмосферы. Приведём несколько замечаний по поводу специфики рассматриваемых задач. Характерные пространственно-временные масштабы процессов различаются на порядки. Размерность современных моделей переноса и трансформации атмосферных примесей достаточно высока и продолжает увеличиваться. При большом числе степеней свободы в моделях данных измерений для нахождения всех неизвестных чаще всего недостаточно, т. е. решение традиционной задачи обращения оператора наблюдений в большинстве случаев неединственно. Поэтому в соответствии с поставленными целями авторы организуют процесс моделирования так, чтобы получать решение в "реальном" времени с использованием доступной на данный момент информации. Такие решения востребованы, например, при ситуациях техногенных и природных катастроф. Целесообразно также предусмотреть возможности использования алгоритмов усвоения как в новых, так и в уже существующих моделях процессов (и в их программных реализациях).
1. Постановка задачи усвоения данных
Рассмотрим прямоугольную пространственно-временную область
X =(Ж1,Ж2,Ж3) е П=[0Л] х [0, /2] X [0,1з] , Ь е [0,Т] , Пт := П х [0,Т] € К4
с границей 8Пт = ^П х [0,Т]. В этой области определим математическую модель типа конвекции-диффузии, которая описывает процессы переноса в атмосфере различных субстанций (примесей, тепла, влаги, радиации и др.):
Ьф = + а1у(иф(х,Ь) - ^(х,^гааф(х,Ь)) = ¡а(х,Ь) + г(х,Ь), (х,Ь) е Пт, (1)
в(х,Ь)ф(х,Ь) = да(х,Ь), (х,Ь) е дПт, (2)
о п
ф(х, 0) = фа(х), х е П, (3)
где ф(х, Ь) — функция состояния; и(х,Ь) = (и^х,Ь),и2(х,Ь),и3(х,Ь)) — вектор скорости ветра; ^(х,Ь) = diag(^1(x, Ь), ^2(х, ¿),^3(х, ¿)) — диагональный тензор коэффициентов диффузии; п — направление внешней нормали к границе области; /а(х, Ь), да(х, ¿), Ф2(х) — априорные значения источников и начальных данных; г(х, ¿) — функция управления (неопределённости), которая вводится в жёсткую структуру модели для усвоения данных. Эту функцию можно рассматривать как объект оценки точности модели в целом и как инструмент уточнения априорно задаваемой интенсивности источников по результатам наблюдений (её смысл будет обсуждаться далее). Задачу определения функций ф из (1)-(3) по известным /а, да, ф°а, г будем называть прямой. Точным решением ф назовём решение прямой задачи, соответствующее некоторому набору функций /, д, ф0, г. Определим оператор наблюдений Н, связывающий функцию состояния модели с результатами наблюдений:
Ф(Ь) = Н(ф(.,¿)) + п(ь), ь е [0,Т] , (4)
где Ф(Ь) е Км — данные наблюдений, М — число измерений в каждый момент времени, П е Км — функция из некоторого множества допустимых значений, представляющая оценку ошибок результатов и моделей наблюдений. Погрешность п в модели измерений предполагается ограниченной во (взвешенной) евклидовой норме в пространстве наблюдений ||п|1м — ^п. Полагаем, что все функции и параметры модели обладают достаточной степенью гладкости, чтобы существовало решение и имели смысл все дальнейшие преобразования. Задачу определения ф(.,Ь*) по (1)-(3) и (4) и заданным при 0 < Ь — Ь* функциям /а, да, фа, Ф назовём задачей усвоения данных.
Такое определение задачи усвоения данных допускает множественные решения. Фактически его дальнейшее уточнение позволяет строить разнообразные подходы к решению задачи. Кратко представим часть из них.
1. Динамико-стохастический подход оценивания параметров методами взвешенных наименьших квадратов, работающий со статистическими моделями [1]:
а) фильтры типа фильтра Калмана [2]; в них функции управления предполагаются гауссовыми случайными величинами, дисперсии которых оцениваются по данным измерений;
б) фильтры частиц [3]. Для различных значений функций управления генерируются множества (ансамбли) решений прямых задач. Решение задачи усвоения выбирается как взвешенная сумма по ансамблю. Веса решения выбираются в соответствии с близостью к данным измерений.
2. Вариационный подход. В этом случае решение задачи усвоения разыскивается как минимум некоторого целевого функционала при условиях, что модели процессов участвуют как ограничения и связи между функциями состояния, параметрами модели и управляющими функциями. Вариационный метод усвоения данных был предложен в 1976 г. [4] (см. также [5]). Универсальность и эффективность этого метода была продемонстрирована на примере четырёхмерной нестационарной модели гидротермодинамики атмосферы для Северного полушария Земли. С тех пор в рамках вариационного подхода разработано много модификаций и новых схем. Можно выделить несколько модификаций вариационных алгоритмов, различающихся тем, как формулируются целевые функционалы и вариационные принципы (со слабыми или сильными ограничениями), в какие части математической модели добавлены функции управления: в начальные данные, в правую часть, в параметры и т.д. [6-8]. Известные модификации ЗВУЛИ и 4БУЛИ различаются интервалами времени, в течение которого учитываются
данные измерений для уточнения решения задачи. При использовании 4DVAR возникает необходимость в итерациях алгоритма в рамках окна усвоения, т. е. данный подход является вычислительно затратным.
3. Представленные выше подходы персекаются с так называемым "nudging" методом — методом усвоения данных, основанным на технике релаксации по Ньютону [9]. В этом методе невязки между измеренными и рассчитанными значениями функции состояния включаются непосредственно в уравнения с некоторыми феноменологическими параметрами релаксации.
Детальные обзоры публикаций (их более 500) с анализом основных подходов к усвоению данных, включая и вариационные методы, представлены в [8, 10-13].
В настоящей работе используется вариационный подход с усвоением на каждом шаге по времени. Рассмотрим некоторую неявную конечномерную аппроксимацию динамической модели (1)—(3):
Ьфj = ф7-1 + т fj, (5)
где j > 1 — номер шага по времени; ф7, ф7 Е RN — дискретизация функции состояния модели на j-м шаге по времени; N — размерность дискретной функции состояния модели; т — шаг сетки по времени; Ь — аппроксимация оператора модели (1)—(3), fj — источники. Для решения задачи усвоения данных часто применяется "явная" схема
ф7 = ф7 + K (иф7 - Ф
где K — некоторая матрица, Я — матричное представление оператора измерений, Ф7 — данные измерений на j-м шаге по времени. В этой схеме функция управления добавляется непосредственно к решению задачи. Данная схема может быть получена также из условий минимума функционала Тихонова (такое представление имеется, например,
в [1, 14]):
2 _ 2 Ф(ф7) = Я(фj - фj) + Яфj - Ф7 + а ф7 - фj ,
M K
где H'llK — норма в некотором весовом пространстве, а — параметр усвоения. Методы такого вида называются 3DVAR. К достоинствам указанной схемы относится то, что для её реализации достаточно иметь результат решения задачи (5) методом прямого моделирования. Сложность в этом случае представляет выбор весовой нормы в регуляризаторе, определяемой в соответствии с пространственной структурой решения. Последнее можно осуществить, например, с применением ансамблевых методов, т. е. запуском большого числа прямых задач для того, чтобы определить, какие части функции состояния при изменении функций управления модели изменяются согласованно. В этом случае веса для нормы рассчитываются как обратные к матрицам ковариации ошибок модели. Вычисление таких матриц сопряжено с рядом сложностей, связанных с необходимостью генерации больших ансамблей для обеспечения репрезентативности. С другой стороны, пренебрежение этими сложностями приводит к тому, что если, например, в качестве нормы рассмотреть сеточную норму Ь2, то в результате шага усвоения модификации будут подвергаться только значения в точках наблюдений, которые не согласуются с физическим смыслом решения: на фоне гладких полей концентраций, порождаемых диффузионными процессами, появляются игольчатые возмущения. Частично проблема гладкости решается добавлением управления не на каждом шаге, а только на некоторых из них (схемы 4DVAR c жёсткими ограничениями). В этом
случае целевой функционал учитывает не один шаг по времени, а несколько (так называемое окно усвоения). Для минимизации такого функционала обычно необходимо проведение итераций.
Рассмотрим другую постановку задачи условной оптимизации, когда на основе данных измерений модифицируется управляющая функция rj в правой части математической модели [15], что затем приводит к модификации решения прямой задачи фр. Такой поход представляется перспективным в задачах усвоения химических данных. В современных обзорах по усвоению данных атмосферной химии (см., например, [8, 13, 16]) отмечается, что основным источником неопределённости в моделях атмосферной химии являются источники выбросов, краевые условия и коэффициенты моделей, которые определяют текущее состояние системы. Начальные данные имеют меньшее значение, так как они под воздействием указанных внешних факторов постепенно "забываются" системой. Вместе с тем исторически широкоизвестные алгоритмы усвоения данных разрабатывались для усвоения метеорологических данных, где количество внешних факторов невелико, и поэтому в них первостепенное внимание уделялось непосредственному уточнению начальных и текущих состояний. В настоящей работе также рассматривается постановка, когда данные измерений используются для уточнения "источников" (emission rate) в модели конвекции-диффузии:
Ьфj = ф-1 + т fj + т rj. (6)
Решением задачи усвоения данных назовём решение задачи минимизации функционала
Ф^, rj) = ||Hфj - ||M + a I|rj ||N
при условии, что выполняется (6). Как и ранее, а — параметр усвоения. Вводя множители Лагранжа ф*3, стационарную точку расширенного функционала
Ф^', rj, ф*3) = ^Hfi - УМ + a ||rj fN + (Ьфj - (ф'-1 + тfj + rrj) , ф*3)
N I yr I I ) > v IN
можно найти из матричной системы
Ь __L \ фj \ фj-1 + тfj
2a
2H *H Ь*
ф*3 ) \ 2H
в которой учтено уравнение, связывающее функции г и ф*3. Здесь Ь* означает оператор, сопряжённый к Ь, <, > — скалярное произведение. Если из данной системы исключить ф*3, то ф3 можно найти из уравнения
(.Н*Н + аЬ*Ь) ф3 = Н*Ф3 + аЬ*(ф3-1 + Л3). В этом случае ф3 является минимумом тихоновского функционала (см., например, [17])
Ф (ф3) = ||Нф3 - Ф3|\М + а ||Ьф3 - (ф3-1 + тР. (7)
Выясним связь между явной и неявной схемами усвоения. Если использовать обозначение (5), то тихоновский функционал (7) приводится к виду
2 _
Ф(ф3)= Н(ф3 - ф3) + Нф3 - Ф3 + а{Ь*Ь(ф3 - ф3), ф3 - ф>)„. (8)
Иными словами, неявная схема эквивалентна явной схеме в весовом пространстве, определяемом оператором Ь*Ь. В дискретном представлении это матрица размерности N х N, где N — число сеточных значений функций состояния модели. Непосредственно вычислять (и хранить) такую матрицу достаточно сложно. Оценивать её с помощью ансамблевых методов также проблематично.
Однако представление (7) позволяет воспользоваться теорией Тихонова для установления свойств решения задачи усвоения в зависимости от значения параметра усвоения а. Для этого обозначим минимум функционала Ф(ф3) через ф3(а). Теорема. Пусть (аналогично [18])
в (а) := ||Нф3 (а) - Ф||М , £(а) := ||Ьф3 (а) - (ф3-1 + тР , 7 (а) := Ф(ф3 (а)) = в (а) + а£(а)
и Н, Ь — матрицы, причём Ь обратима.
— Если увеличивается параметр усвоения, то решение приближается к решению задачи прямого моделирования:
I 2
I М '
lim ¿(а) = 0, lim J(а) = lim ß(а) = ||ЯХ-1(ф'-1 + тfj) - Ф7
а^+те а^+те а^+те
— Если параметр усвоения уменьшается, то данные измерений (и шум в них) "воспроизводятся" точнее в решении задачи усвоения:
lim а£ (а) = 0, lim J (а) = lim в (а) = 0.
— Функции в (а), ¿(а), J (а) непрерывны при а > 0.
— Функция J(а) выпукла и дифференцируема J'^) = ¿(а).
— Функция ¿ (а) монотонно невозрастающая, а функции J (а), в (а) монотонно неубывающие при а > 0.
Доказательство проводится по схеме доказательства леммы 3 из работы [18, с. 22-24], где функционал ¿(а) взят в более простом виде ¿(а) := ||ф7(а)||м. Результат является естественным обобщением данной леммы.
Фактически параметр усвоения регулирует невязку между данными измерений и результатом усвоения данных. Если данные неточны, эта невязка не обязательно должна быть нулевой. Для выбора параметра усвоения (регуляризации) в теории некорректных задач применяется принцип невязки Морозова, в соответствии с которым параметр следует выбирать так, чтобы невязка измерений в (а) соответствовала погрешности измерений.
Чтобы удовлетворить условию о величине погрешности измерений, можно воспользоваться монотонностью функционала и предложить следующий алгоритм для выбора параметра усвоения (регуляризации):
— выберем целевой уровень невязки 8*
0 < ||Hфj - Ф7||M « 8* < ||Я£-1(ф7-1 + Tfj) - Ф^||м ;
— параметр усвоения вычисляется итерационным алгоритмом
8*
а1 = 1, ак+1 := м —-- || ак.
Для выбора способа оценки 8* уточним вид оператора наблюдений. Предположим, что в момент времени t проводятся измерения в M точках области {Xm Е П|ш = 1, M} точного решения ф и данные измерений содержат нормальный шум с известными стандартными отклонениями {om > 0|m = 1,M}. Тогда оператор измерения имеет вид Н(ф) = {ф(х1, t) + tfi£i(t),... ,ф(Хм (t),t) + ом СМ (t)}*, где * означает операцию транспонирования, а £m — нормально распределённые случайные величины из N(0,1). Без ограничения общности предположим, что в m-й точке на j-м шаге измеряются значения вектор-функции состояния модели с индексом i3m. Определим весовую норму
в пространстве измерений
1Ф1
м
\
м
£
m=1
—— ] . Если данные точны (om = 0), то
Om
в качестве нормы ||.||м выберем евклидову норму. Тогда для оценки 8* можно использовать следующую лемму.
Лемма. Если данные измерений — это точное решение в некотором наборе точек, возмущённое гауссовым шумом
" N (0,1), т = Т~М,
—j
m
ф .j + оmСm, Сп
то с вероятностью p
8
\
м
£
m=1
'—m - ф.
<
о„
Xinv
(M,p),
где Хл™(М,р) — решение уравнения Р(хМ < Хл™(М,р)) = р, Хм — случайная величина с распределением хи-квадрат с М степенями свободы.
Доказательство. По определению результат измерения и случайная величина с распределением хи-квадрат связаны соотношением
м
£
m=1
—m ф
м
Om
Ее.
m=1
m = Хм,
(9)
откуда и получаем утверждение леммы. □
Фактически для параметра усвоения вводится новая шкала: в качестве оценки ошибки "выбирается" вероятность р её оценки "сверху". Чем больше эта вероятность, тем более завышенной (надёжной) становится оценка.
Добавление оператора прямой задачи в регуляризатор уменьшает "шоковые" эффекты от включения данных в модель (сглаживает решение). Недостатком данной схемы является необходимость решать систему прямой и сопряжённой задач, сгенерировать которую для "готовых" моделей может быть сложно. В общем виде решение системы требует привлечения итерационных методов, что снижает производительность алгоритма. Цель настоящей работы — представить вычислительно эффективную технологию, решающую проблему усвоения данных.
2
2
m
2. Вычислительная технология мелкозернистого вариационного усвоения данных
Для решения многомерной прямой задачи применяется метод расщепления (по пространственным переменным) [19]. Воспользуемся этим методом и для решения задачи усвоения данных [5, 15]. Чтобы задействовать технологии массового параллелизма,
будем использовать аддитивно-усреднённые схемы расщепления (аналогично [20]) и, следуя [21], схемы декомпозиции. Введём на временном отрезке сетку по времени:
^ = {0 < ... < ... = Т}.
На каждом из интервалов времени общую модель конвекции-диффузии (1)-(3) в областях П х [¿•7-1,£7], ] > 1 можно представить в виде схемы расщепления, имеющей свойства суммарной аппроксимации, этапы которой представляют локально одномерные модели процессов конвекции-диффузии
вл^ + д(в>мж(х.^)) _ а ( (х )аф^ = Л(х,г) + Гк(х,г), (10)
от ахк ахк
(х. ¿) е П х (V-1 ] . (11)
^ (х. *) + в (х. 4)фк (х. *) = 0(х, *). (12)
дп
(х.*) е (<Ш П ({Хк = 0} и {Хк = /к})) х р-1.^'] . (13)
фк (х.*7-1) = ф (х.*7-1) . ^ > 1. 7к > 0. к =1. 2. 3. (14)
3 3
где £ 7к = 1,5^ /к = /. Приближённое решение исходной задачи в момент времени к=1 к=1 * = вычисляется как
3
ф(х.*7) = ^7кфк(х.*7). х е П. (15)
к=1
При заданном параметре усвоения для каждой одномерной задачи в схеме расщепления многомерной модели конвекции-диффузии можно построить прямой вычислительно эффективный алгоритм, не требующий итераций.
2.1. Неявный алгоритм усвоения данных для одномерной модели конвекции-диффузии
Каждая нестационарная одномерная подмодель конвекции-диффузии из (10)-(14) аппроксимируется на пространственно-временной сетке так, чтобы получить трёхточеч-ную систему линейных алгебраических уравнений
_а*$+1 + 6гф7 _ сгф7-1 = ф'-1 + т(/ + г)7'. г = 0^. со = 0. = 0. (16)
где верхний индекс ] > 1 означает номер шага по времени, нижний г — номер узла пространственной сетки; N — число интервалов пространственной сетки; т = А*7 — шаг по времени; ai.bi.Ci — коэффициенты численной схемы. Для краткости номер-индекс к этапа схемы расщепления опустим. Трёхточечные аппроксимации можно получить, например, с применением конечно-разностной схемы направленных разностей или дискретно-аналитической схемы [22, 23].
Запишем с помощью матричного оператора Ь систему (16) в виде
Ьф7 = ф'-1 + т (Р + г7).
Для удобства введём обозначения: множество индексов узлов сетки на ^'-м шаге, в которых производятся измерения = {г^! т =1. М}, индикаторную функцию системы
измерений, функцию продолжения нулём данных измерений из точек измерений на всю область, функцию продолжения единицей значения стандартных отклонений соответствующих измерений:
0, г ^ I3,
Фт, Зт|г гт 0, г ^ ,
с г
ст, Зт|г = г^ 1, г / Р.
:17)
Решение задачи получим из условия минимума расширенного функционала, принимая в качестве ограничений уравнения (16) схемы решения прямой задачи. Вводя множители Лагранжа (сопряжённые функции), определим расширенный функционал со слабыми ограничениями в виде
N
1 / " / _ Ф(ф3,г3,ф*3) = 1(Х] ' Фг Ф 12 т3
N
2
, г=0
С
+«Е
ГЛ 2
г=0
N
+ К (-аФЗ+1 + Ьф3 - сгф?-! - Ф3-1 - т(Ц + г3)) ф.
г=0
%.
:18)
Записывая условия стационарности функционала (18) к вариациям функций ф3, г3, ф*3, получим:
a) систему основных уравнений дф3 Ф(ф3, г3 , ф*3) = 0, эквивалентную схеме для прямой задачи (16);
b) уравнения для оценки неопределённостей 8гз Ф(ф3, г3, ф*3) = 0 в виде
аг3 - ф*3 = 0, г = 0, N;
19)
с) систему сопряжённых уравнений дф Ф(ф3, г3, ф*3) = 0, эквивалентную следующей системе трёхточечных уравнений:
-С+1ф*3+1 + Ьф*3 - аг-1ф*3-1 = -13 (ф - Ф') т/(ас2),
г = 0, N, а-1 = с.м+1 = 0.
(20)
Используя уравнения (19), системы прямой и сопряжённой задач объединим в систему матричных уравнений [15, 24]
-АгФ3+1 + ВФ3 - СФ3-1 = ^, г = 0, N
где
А
а 0
0 Сг+1 0,
г = - 1,
г = N
С
В
Ьг
ьг -Т
т/(аСт2) Ьг
Ф
3
ф3 ф*3
к3
0,
С 0
0 <а -1
фТ1
т /(ас 2)
1, N
г = 0, N,
которая решается методом матричной прогонки.
Пример решения одномерной задачи усвоения данных представлен на рис. 1. Условия эксперимента: решается задача переноса-диффузии с заданными начальными данными при £ = 0 без источника. Решение этой задачи принимается в качестве "точного".
0
а б в
Рис. 1. Результат усвоения данных для различных направлений "ветра". Верхний ряд — точное решение, нижний — результат усвоения; у = 0.025, и = -0.5 (а), и = 0 (б), и = 0.5 (в)
Далее задача усвоения решается с нулевыми начальными условиями также без источника. В роли "наблюдений" выступают наборы значений точного решения в некоторых фиксированных точках области П (стационарная сеть). Таким образом, пространственно-временная структура решения восстанавливается только по данным наблюдений. Графики (см. рис. 1) иллюстрируют вид "точного" решения при трёх вариантах направления скорости переноса (верхний ряд) и решения задач усвоения для соответствующих условий (нижний ряд). Из рисунков видно, что структура решений восстанавливается физически адекватно, но при этом проявляются эффекты, обусловленные "точечной" структурой наблюдений. Отметим: чтобы получить аналогичные решения методом ЗВУАИ,, пришлось бы в соответствии с (8) дополнительно вычислять матри-
N
цу ковариации. Если кроме слагаемого целевого функционала (г7), отвечающего
г=0
за малость управлений, добавить ещё слагаемое, обеспечивающее малость его сеточных производных (т.е. вместо сеточной нормы Ь2 использовать, например, сеточную норму Н1), то гладкость решения можно ещё более увеличить.
2.2. Мелкозернистое усвоение в многомерном случае
Поясним термин "мелкозернистое усвоение" через различие подхода к комбинированию расщепления и усвоения данных, развитого в работах [5, 25, 26], и подхода, разрабатываемого в [4, 15, 24] и в настоящей статье. В первом случае для отдельных стадий расщепления по физическим и пространственным переменным строятся отдельные сопряжённые задачи, но цикл усвоения проводится по отношению ко всему шагу по времени. Сопряжённые задачи для шагов расщепления оказываются связанными, и для решения системы оптимальности для целевого функционала необходимо проведение
итерации. Во втором подходе для отдельных шагов строятся не только сопряжённые задачи, но и сами целевые функционалы и соответствующие им системы оптимальности. При этом сопряжённые задачи для различных этапов расщепления оказываются несвязанными. Если в этом случае найти эффективный метод решения частных задач усвоения, то можно получить высокоэффективный в вычислительном плане алгоритм усвоения данных.
Введём в области П сетку
3
П = До, о = {0 < ... < х < ... < = 4}
с числом узлов N в к-м измерении. Узел в такой трёхмерной сетке будет определяться мультииндексом 2 = {¿1,22,23}. Для каждого этапа расщепления будем рассматривать множества узлов, соответствующих линиям сетки в заданном измерении. Таким образом, линия сетки определяется значениями мультииндексов на всех позициях за исключением к-й. Значение к-го мультииндекса будет определять элемент из данной линии сетки. Далее на линиях сетки, свободных от точек наблюдения, будем решать одномерную прямую задачу, соответствующую этапу расщепления, а на линиях сетки, где располагаются точки измерений, прямой вариант модели заменим на задачу нахождения стационарных точек соответствующих расширенных функционалов. Таким образом, каждое измерение оказывается на пересечении линий сетки по числу пространственных измерений и усваивается соответствующее количество раз. Аналогично одномерному случаю (17) продолжим значения измерений на всю область: введём функцию-индикатор системы измерений , продолжение нулём поля наблюдений Ф-7 и продолжение единицей значения стандартных отклонений ошибок измерений а:
Фк(Фк, г-, Фк) = 2 + ак £ (7)2^ т+
+ Фк - Ф7-1 - т(/ - г-), Фк>к , з > 1, к = М, (21)
где — дискретизация оператора расщеплённой прямой задачи (10)-(14), имеющая вид, аналогичный левой части (16), (.,.)& — скалярное произведение в пространствах . Отметим, что параметр усвоения выбирается для каждой координатной линии. Все подзадачи имеют параллельную структуру и могут решаться одновременно. После их решения по алгоритму (см. раздел 2.1) окончательно результат задачи усвоения данных на момент времени £ = получаем по алгоритму (15), т.е. результирующие поля Фк суммируются с весами .
3. Численное исследование устойчивости
Важнейшим вопросом, возникающим при решении обратных задач, является вопрос чувствительности и устойчивости решения к ошибкам в данных измерений. Из теоремы (см. с. 74) можно заключить, что эта устойчивость зависит и от параметра усвоения а.
На примере двухмерной версии модели (1)-(4) и соответствующей ей схемы (10)—(15) рассмотрим, как параметр усвоения влияет на точность решения при нескольких вариантах задания погрешности в данных наблюдений. Численные эксперименты
по усвоению организованы так же, как и в разделе 2.1, т.е. без начальных данных и с нулевыми источниками. В таблице представлено описание параметров системы измерений. Параметры модели: и = {0.5, 0.5}, ^ = 0.001, сеток: N = 100, Т =1, Ых = 100, 1х = 1, Ыу = 100, 1у = 1. На рис. 2 приведено сравнение решения задачи усвоения данных (а) с модельным "точным" решением (б) в заданный момент времени.
Параметры системы измерений при проведении численных экспериментов*
Параметр т
системы измерений 1 2 3 4 5 6 7 8 9 10 11 12
^Хт 33 33 67 67 25 25 75 75 40 60 40 60
ЪУт 33 67 33 67 25 75 25 75 60 40 40 60
от (слева) 0.1 1. 0.5 1. 1. 2. 1. 0.5 1. 0.5 3. 0.1
от (справа) 4. 40. 20. 40. 40. 80. 40. 20. 40. 20. 120. 4.
*т — номер поста наблюдения, гХт — индекс узла сетки по оси X, где происходят измерения, гУт — индекс узла сетки по оси У, где происходят измерения, ат — дисперсия соответствующего измерения
Рис. 2. Решение задачи усвоения данных в заданный момент времени (Ь = 15) с двенадцатью регулярно расположенными измерительными инструментами (а), "точное решение" (б)
Рис. 3. Зависимость ошибки решения от параметра усвоения р при "малых" (а) и "больших" (б) значениях дисперсии ошибки измерений (см. таблицу)
Рис. 4. Зависимость ошибки решения от числа постов наблюдений М при "малых" (а) и "больших" (б) значениях дисперсии ошибки измерений (см. таблицу)
а
В первом численном эксперименте (рис. 3) изучается зависимость абсолютной ошибки решения задачи усвоения данных от параметра усвоения р. Для данных с "малыми" ошибками их учёт приводит к улучшению результатов, а в случае данных с шумом лучший результат получается при средних значениях параметра усвоения.
На рис. 4 представлена зависимость абсолютной ошибки решения задачи усвоения данных от числа постов наблюдений. Параметр усвоения р = 0.1. В каждом эксперименте добавляются по четыре дополнительных поста М. Из приведённых графиков можно заключить, что при фиксированном значении параметра усвоения р в случае точных данных увеличение числа постов наблюдений приводит к улучшению решения, а в случае данных с шумом дополнительные посты приводят к ухудшению решения.
Заключение
В работе представлен алгоритм усвоения данных на одном шаге по времени в дискретном аналоге динамической модели процессов. В отличие от стандартных ("явных") схем ЗВУЛИ в предложенном алгоритме нет разделения на вычисление "первого приближения" и последующее вычисление матриц ковариации ошибок первого приближения. Фактически эти матрицы указывают, какие части функции состояния изменяются согласованно. В случае моделей большой размерности вычисление и хранение таких матриц с размерностью, равной квадрату числа переменных, представляет определённую проблему. Вместо этого через оптимизацию целевого функционала находится не сама функция состояния, а некоторая управляющая функция, которая добавляется в основную модель для обеспечения необходимых для усвоения данных степеней свободы. Решение задачи усвоения данных получается в результате прямого шага с учётом найденных значений управляющей функции. Тем самым происходит автоматическое согласование изменений различных частей функции состояния. В работе показана связь между предложенным подходом и схемой ЗВУЛИ.
При решении многомерных задач вариационное усвоение данных производится на базе схем расщепления и приводит к вычислительно эффективным прямым алгоритмам для уравнений конвекции-диффузии в случае усвоения контактных измерений. Вычислительная эффективность алгоритма обусловлена использованием схемы усвоения на отдельных этапах расщепления и наличием эффективной процедуры решения
системы оптимальности для каждого элементарного шага расщепления. Применение вариационных принципов со слабыми ограничениями позволяет перенести работу по усвоению данных и управлению неопределённостями на уровень схем декомпозиции и расщепления функционалов и моделей.
Эффективность предложенного алгоритма подтверждается численными экспериментами. Более строгое теоретическое обоснование представленного материала может быть предметом будущей работы.
Список литературы
[1] Брдйсон А., Ю-Ши Хо. Прикладная теория оптимального управления. М.: Мир, 1972. 544 с.
[2] Kalman R.E. A new approach to linear filtering and prediction problems // Trans. of the ASME-J. of Basic Eng. 1960. Vol. 82, ser. D. P. 35-45.
[3] Gordon N.J., Salmond D.J., Smith A.F.M. Novel approach to nonlinear/non-Gaussian Bayesian state estimation // IEE-Proceedings-F. 1993. Vol. 140. P. 107-113.
[4] Пененко В.В., Образцов Н.Н. Вариационный метод согласования полей метеорологических элементов // Метеорология и гидрология. 1976. № 11. C. 1-11.
[5] Пененко В.В. Методы численного моделирования атмосферных процессов. Л.: Ги-дрометеоиздат, 1981. 352 с.
[6] Dimet F.-X. L., Talagrand O. Variational algorithms for analysis and assimilation of meteorological observations: Theoretical aspects // Tellus. A. 1986. Vol. 38A. P. 97-110.
[7] Talagrand O., Courtier P. Variational assimilation of meteorological observations with the adjoint vorticity equation // Quart. J. Roy. Meteor. Soc. 1986. Vol. 113, No. I: Theory. P. 1311-1328.
[8] Elbern H., Strunk A., Schmidt H., Talagrand O. Emission rate and chemical state estimation by 4-dimensional variational inversion // Atmos. Chem. Phys. 2007. Vol. 7. P. 3749-3769.
[9] Anthes R.A. Data assimilation and initialization of hurricane prediction models //J. Atmos. Sci. 1974. Vol. 31. P. 702-719.
[10] Blum J., Le Dimet F.-X., Navon I.M. Data assimilation for geophysical fluids // Handbook of Numerical Analysis. Els. Sci. 2008. Vol. 14: Special Volume. P. 377-434.
[11] Navon I.M. Data assimilation for numerical weather prediction: A review // Data Assimilation for Atmospheric, Oceanic and Hydrologic Appl. 2009. Vol. 1. P. 21-65.
[12] Data Assimilation: Making Sense of Observations / Eds. W. Lahoz, B. Khttatov, R. Menard. Heildelberg: Springer, 2010. 718 p.
[13] Sandu A., Tianfeng C. Chemical data assimilation — an overview // Atmosphere. 2011. Vol. 2. P. 426-463.
[14] Freitag M.A., Potthast R.W.E. Synergy of inverse problems and data assimilation techniques // Large Scale Inverse Problems Comput. Methods and Appl. in the Earth Sci.: Radon Series on Comput. and Appl. Math. 2013. Vol. 13. P. 1-54.
[15] Пененко В.В. Вариационные методы усвоения данных и обратные задачи для изучения атмосферы, океана и окружающей среды // Сиб. журн. вычисл. математики. 2009. Т. 12, № 4. C. 341-351.
[16] Baklanov A., SchlUnzen K., Suppan P., Baldasano J. et al. Online coupled regional meteorology chemistry models in Europe: Current status and prospects // Atmos. Chem. Phys. 2014. Vol. 14. P. 317-398.
[17] Агошков В.И. Теория и методы решения задач вариационной ассимиляции образов. M.: ИВМ РАН, 2012. 147 с.
[18] Тихонов А.Н., Гончлрский А.В., Степанов В.В., ЯголА А.Г. Численные методы решения некорректных задач. M.: Наука, 1990. 230 с.
[19] Марчук Г.И. Методы вычислительной математики. M.: Наука, 1989. 608 с.
[20] Самарский А.А., Вабищевич П.Н. Аддитивные схемы для задач математической физики. M.: Наука, 2001. 306 с.
[21] Бенсусан А., Лионе Ж.-Л., Темам Р. Методы декомпозиции, децентрализации, координации и их приложения // Методы вычислительной математики. Новосибирск: Наука, 1975. С. 144-270.
[22] Penenko V.V., Tsvetova E.A. Discrete-analytical methods for the implementation of variational principles in environmental applications // J. of Comput. and Appl. Math. 2009. Vol. 226. P. 319-330.
[23] Пененко А.В. Дискретно-аналитические схемы для решения обратной коэффициентной задачи теплопроводности слоистых сред градиентными методами // Сиб. журн. вычисл. математики. 2012. Т. 15, № 4. C. 393-408.
[24] Пененко А.В. Некоторые теоретические и прикладные вопросы последовательного вариационного усвоения данных // Вычисл. технологии. 2006. Т. 11. Спецвыпуск: Избранные докл. Междунар. конф. по измерениям, моделированию и информационным системам для изучения окружающей среды. ENVIR0MIS-2006. Томск, 2006. Ч. 2. C. 35-40.
[25] Marchuk G.I., Zalesny V.B. A numerical technique for geophysical data assimilation problems using Pontryagin's principle and splitting-up method // Rus. J. of Numerical Analysis and Math. Modelling. 1993. Vol. 8, No. 4. C. 311-326.
[26] Агошков В.И., Ипатова В.М., Залесный В.Б. и др. Задачи вариационной ассимиляции данных наблюдений для моделей общей циркуляции океана и методы их решения // Изв. РАН. Физика атмосферы и океана. 2010. Т. 46, № 6. C. 734-770.
Поступила в 'редакцию 22 февраля 2014 г., с доработки — 29 мая 2014 г.
