Применени ансамблевого пи-алгоритма в задаче обратного моделирования Текст научной статьи по специальности «Математика»

АПВПМ-2019

ПРИМЕНЕНИ АНСАМБЛЕВОГО ПИ-АЛГОРИТМА В ЗАДАЧЕ ОБРАТНОГО МОДЕЛИРОВАНИЯ

Е, Г, Климова

Институт вычислительных технологий СО РАН, 630090, Новосибирск

УДК 551.509.313

DOI: 10.24411/9999-016А-2019-10038

Решение задачи обратного моделирования для заданной модели процесса и набора данных наблюдений циклически по времени является задачей усвоения данных. В проблему обратного моделирования включена также задача оценки параметров модели. Использование динамико-стохастического подхода к решению задачи усвоения данных привело к развитию ансамблевых методов: ансамблевого фильтра Калмана, ансамблевого сглаживания Калмана. В задаче усвоения данных при моделировании процессов в окружающей среде оцениваемый вектор и вектор данных наблюдений имеют высокую размерность, что требует разработки эффективных алгоритмов. Ансамблевый пи-алгоритм разработан как эффективная реализация стохастического ансамблевого фильтра Калмана. Он может быть использован в задачах оценки состояния окружающей среды, а также параметров модели. Одной из проблем ансамблевых фильтров является возможная расходимость алгоритма со временем. В докладе рассматриваются подходы к практической реализации ансамблевого фильтра Калмана и ансамблевого сглаживания Калмана с использованием пи-алгоритма, связанные с адаптивной коррекцией поведения алгоритма по времени. Приводятся результаты численных экспериментов с 1-мерной тестовой моделью. Ключевые слова: усвоение данных, ансамблевые алроритмы.

Введение

Решение задачи обратного моделирования для заданной модели процесса и набора данных наблюдений циклически по времени является задачей усвоения данных. В проблему обратного моделирования включена также задача оценки параметров модели. Для ее решения часто применяется байесовский подход, при этом данные наблюдений и оцениваемые параметры считаются случайными величинами. В случае, если рассматриваемые случайные величины являются гауссовскими, а модели прогноза и наблюдений линейными, такая постановка задачи эквивалентна вариационной постановке задачи усвоения данных (4DVAR) [7]. Для аппроксимации рассматриваемых в алгоритме ковариационных матриц может быть использован ансамбль (выборка) прогнозов и наблюдений [4]. В случае, если оптимальная оценка ищется в конце заданного временного интервала, то задача сводится к ансамблевому фильтру Калмана. Если рассматривается оптимальная оценка на заданном временном интервале, такая задача называется задачей ансамблевого сглаживания [4]. В задаче усвоения данных при моделировании процессов в окружающей среде оцениваемый вектор и вектор данных наблюдений имеют высокую размерность, что требует разработки эффективных алгоритмов.

Одной из проблем ансамблевых фильтров является возможная расходимость алгоритма со временем. В обзоре [8] перечисляются основные подходы к улучшению сходимости ансамблевого фильтра Калмана. Одним из популярных способов регулирования сходимости ансамблевого фильтра Калмана является использование так называемого "inflation factor" (увеличивающий множитель). В работе [8] отмечаются такие способы борьбы с расходимостью ансамблевого фильтра Калмана, как увеличивающий множитель (multiplicative inflation) и увеличивающее дополнительное возмущение (additive inflation). В нелинейном случае в настоящее время применяются итерационные алгоритмы сглаживания, представляющие собой вариационную задачу с применением ансамблей для вычисления ковариаций и линеаризации операторов прогноза и наблюдений [3,5]. В строго нелинейном негауссовском случае используется метод частиц, который основан на байесовском подходе [3]. В данных подходах также рассматриваются способы управления сходимостью алгоритмов. Ансамблевый пи-алгоритм эффективен при использовании при реализации стохастического ансамблевого

ISBN 978-5-901548-42-4

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

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

1 Задача обратного моделирования

Задачей обратного моделирования называют оценку неизвестного вектора параметров большой размерности х € RK по вектору данных наблюдений у G RL [13]. Если предположить, что

у = М (х) + S,

где М — известный оператор прогноза-наблюдений, 5 — случайный шум с заданной функцией распределения, то оценка вектора состояния х то данным наблюдений у производится на основе байесовского подхода. Задачу последовательного оценивания по времени неизвестного величины по временному ряду данных наблюдений называют задачей усвоения данных [13].

Предположим, что изменение по времени оцениваемой величины х^ описывается моделью

Хк+1 = fk+1,k (хк) + Ilk, где к — помер шага по времени. Кроме того, известны данные наблюдений ук'.

Ук = hk(xk) + £к-

В этих формулах щ, Ек — случайные ошибки прогноза и наблюдений, соответственно, fk+1,k — оператор модели, hk — оператор трансформации прогнозируемой переменной в наблюдаемую.

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

p(y]ix)p(x)

РШ = ~МГ

Известны различные варианты оценки состояния по данным и прогнозу: Р(%1 Ьм)> к>1 — прогноз, p(xklVk,i) — фильтрация, р(хк,о1'Ук,1) — сглаживание,

где Xkfi = [хк,Хк-1у • • • , хо], Ук,1 = {Ук, • • • ,у1}- Обозначения и формулировки взяты из обзора [3]. В линейном гауссовском случае решением задачи фильтрации является фильтр Калмана, решением задачи сглаживания — сглаживание Калмана. В работе [4] предложено использовать метод Монте-Карло для решения задач фильтрации и сглаживания, так называемые ансамблевые алгоритмы фильтрации и сглаживания. В ансамблевом фильтре Калмана в нелинейном случае прогноз ковариационной матрицы осуществляется с помощью нелинейной модели, при этом нарушается условие гауссовости. Также в этом случае оценка на шаге анализа является приближением оценки минимальной дисперсии (linear variance minimizing). Для решения задачи сглаживания в нелинейном случае в настоящее время рассматриваются итерационные алгоритмы сглаживания с применением ансамблей. Как отмечается в [5], итерационные методы дают хорошее приближение для слабо нелинейных моделей. В нелинейном негауссовском случае используется метод частиц, который также основан на байесовсом подходе.

2 Применение ансамблевого подхода к задаче оптимальной фильтрации

Стохастический ансамблевый фильтр Калмана [4] состоит из ансамбля прогнозов {х['п, п = 1, • • • , N}

4'" = f («) + ^-1 (1)

и ансамбля анализов п = 1, • • • , N}

хак'п = + Кк[Уп + еп — h(xfk'n)]. (2)

В этих формулах N — количество элементов ансамбля. Ансамбли (1) и (2) задают выборку значений «истины», при этом среднее по выборке значение будет являться оптимальной оценкой, а отклонение от среднего — ансамблем ошибок анализа и прогноза, соответственно. Для осуществления ансамблевого варианта алгоритма фильтра Калмана требуется задание ансамбля ошибок наблюдений [е", п =1, • • • , N} а также ансамбля

ошибок прогноза [dxkn = хкп — х['п, п = 1, • • • , N}, где х['п = N J^N=1 %к'п и ансамбля шумов модели

[v'n-1, п =1, • • • , N} Е['Пп-1 {'П'гп-1)Т] = Qk-1 • Матрица Кк имеет вид

Кк = р[ Щ (Нк р[ Hi + Rk )-1, где Р[ и Rk — матрицы, оцениваемые по ансамблю

1 N т 1 N

pl = jhï^^Ы'п) ,Rk = N^îT.сп(^п)т,

п=1 п=1

Hk ^^^^^^^^^^адный оператор h(х^) относительно xJk'la:

h(xkn) - h(xkn) + Hk dxkn-

Детерминированный ансамблевый фильтр Калмана (шаг анализа) состоит из уравнения для среднего значения

Зг = + Kk [у n — h(xk'n)]

и оценки ансамбля ошибок анализа такой, чтобы соответствующая ковариационная матрица удовлетворяла уравнению фильтра Калмана = (I — KkHk)Р£ [8].

Одним из преимуществ стохастического фильтра Калмана является то, что отклонение элемента ансамбля от среднего (ensemble spread) имеет вид, аналогичный ошибке оценки — отклонению среднего значения от "истины" (ensemble skill) [8,12].

3 Проблемы реализации ансамблевых алгоритмов

Усвоение данных для современных задач моделирования процессов в окружающей среде имеет большую размерность как прогнозируемого вектора, так и вектора наблюдений. Всвязи с этим актуальной является разработка эффективных алгоритмов реализации процедуры усвоения данных. Одним из способов создания экономичного алгоритма является использование методов трансформации ансамбля возмущений. Ансамблевый пи-алгоритм является эффективной реализацией стохастического ансамблевого фильтра Калмана [1,10]. Ансамблевый пи-алгоритм является стохастическим фильтром, в котором ансамбль ошибок анализа получается с помощью матрицы трансформаций, не зависящей от узла сетки. Это делает возможным реализацию алгоритма локально, при этом для реализации ансамблевого пи-алгоритма требуются операции с матрицами порядка размерности ансамбля. Как видно из формул для оценки параметра в ансамблевом алгоритме, матрица трансформаций одна и та же для оценки моделируемой переменной и параметра модели. Обоснование метода и детали его реализации содержатся в [1,10].

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

Малый размер выборки также является одной из причин расходимости алгоритма. Расодимость бывает двух типов. В первом случае происходит рост среднеквадратиеской ошибки. Во втором случае быстро убывает след матрицы ковариаций ошибок оценивания и в результате весовая матрица при данных наблюдений

становится близкой к нулевой. Для предотвращения расходимости второго типа в алгоритмах ансамблевого фильтра Калмана (EnKF), ансамблевого сглаживания Калмана (EnKS) используется «увеличивающий множитель» (inflation factor), при этом либо элементы ковариационной матрицы умножаются на некоторый множитель (multiplicative inflation), либо к ансамблю возмущений добавляется дополнительное возмущение (additive inflation). Следует отметить, что эти приемы используются и в ряде итерационных алгоритмов сглаживания с применением ансамблей, а также в методе частиц с применением ансамблей [8].

Следует отметить, что применение мультипликативного множителя в линейном случае эквивалентно увеличению весового множителя при более близко расположенных к заданному моменту времени данных [9].

В целях улучшения сходимости алгоритмов по времени популярным является адаптивное уточнение следа матрицы ковариаций ошибок оценивания. В основе этого подхода лежит свойство вектора невязок

гЩ:

г = У0 — H(хъ) « £0 - H(£Ь),

где уо — вектор данных наблюдений, хъ — вектор прогностических переменных (первого приближения), H — оператор наблюдений, е0 , еЬ — векторы случайных ошибок прогноза и наблюдений, соответственно. Вектор невязок обладает свойством

ггт = 100Щ + НёЬе[Нт = R + НРЬНТ,

где R, Ръ — матрицы ковариаций ошибок наблюдений и первого приближения, соответственно. Для оценки увеличивающего множителя используется формула

гт г — tr(R) Р = tr(HPbHT) '

при этом принято предполагать, что rTг « rTr, tr означает след матрицы.

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

4 Пример управления сходимостью ансамблевого алгоритма с помощью уеличивающего множителя

Рассмотрим подход к управлению сходимостью ансамблевого фильтра Калмана на примере ансамблевого пи-алгоритма для 1-мерной модели переноса и диффузии пассивной примеси с неизвестным параметром (эмиссией примеси). Для решения уравнения переноса и диффузии применялся полу-лагранжев метод, при этом бралась неявная схема по времени и схема центральных разностей по пространству. Для решения конечно-разностного аналога уравнения диффузии использовался метод циклической прогонки. Уравнение решалось на отрезке (0,1), при этом рассматривались периодические граничные условия. Задавалось 240 узлов сетки, и = 1, к2 =0, 6 х 10-3.

Рассмотрим конечно-разностный аналог уравнения переноса и диффузии в виде

Фк+1 = Ак фк + дк,

где Ак — линейный оператор, к — номер шага по времени. Заданные начальные значения ф0, д0 считались «истинными». Для получения начальных данных для прогноза по модели к «истинным» начальным данным добавлялось возмущение = фг0 + 6, 6 = N(0, So). Через N(a,b) обозначена случайная величина, распределенная по нормальному закону с математическим ожиданием, равным а и дисперсией, равной Ь.

Для организации численных экспериментов задавались: ансамбль начальных полей = + ân, ôn ~ N(0, S0), n=1, •• •, NenS; дЦ = + 5ng, 5ng ~ N(0, dg0), n = 1, •• •, NenS, наблюдения % = Фь(0) + ¿0; ¿0 ~ N(0,^); ансамбль наблюдений с возмущениями уП = у0 + ^П, ^о ~ N(0, е0), n = 1, • • •, Nens. Через Nens обозначено число элементов ансамбля. Данные наблюдений считались известными во всей области интегрирования. Прогноз осуществлялся в течение ^=240 шагов по времени, усвоение данных проводилось па каждом шаге по времени, при этом производилась оценка значений ф и д. Численные эксперименты проводились для значений s0 = е0 = 0.01 Nens = 20 и 40.

Во всех численных экспериментах рассматривался вариант R = ¿01. При анализе в узле сетки I брались данные наблюдений из интервала (I — id,l + id). При этом при анализе в узле сетки I вместо матрицы R

бралась матрица R = R о e-0-5(pil/bc)2, где рц — расстояние между узлом сетки и наблюдением, " о " — знак поэлементного умножения. В проводимых экспериментах брались значения id =5, be = 5Дх (Дх — шаг сетки). Такой алгоритм носит название локализация. Его принято использовать в ансамблевых методах для подавления ложных ковариаций на больших расстояниях, возникающих из-за малого размера выборки (ансамбля).

В численных экспериментах "истинное" значение параметра задавалось кусочно-постоянным в виде дискретного аналога функции д(х, t), где:

д0(х) , 0 <t < 120 0.8 * д0(х), t> 120 '

где

0.1 , 0.375 < х < 0.625 0, 0.375 > х Л х > 0.625 .

Для предотвращения расходимости алгоритма со временем использовался так называемый «увеличивающий множитель» infl (inflationfactor) [8]. То есть, на шаге прогноза ансамбль возмущений домножался на константу: вместо Pf рассматривалась матрица, все элементы которой умножены на infl.

В численных экспериментах использовалась ансамблевый пи-алгоритм стохастический ансамблевый

n

n f = 1 n f = 1 . 3 n

дом узле сетки отдельно, что возможно при реализации пи-алгоритма, поскольку он является локальным

n

рассчитывается как среднее значение по подобласти узла сетки, для которого производится усвоение. В

n

n

inf 11 = ft * infl + (1 — Р) * inf 10,

где inflo — значение с предыдущего шага по времени, Р = 0.8. Результаты экспериментов приведены на рис. 1 и рис. 2. На этих рисунках представлено поведение среднеквадратической ошибки переменных ф и д,

n

д(х, ) =

0(х)

time step

Рис. 1: Среднеквадратическая ошибка оценки прогнозируемой переменной

time step

Рис. 2: Сродноквадратичоская ошибка оценки параметра

Как видно из рисунков, расходимость алгоритма усвоения для заданной модели происходит после изменения "истинного" значения неизвестного параметра. Использование параметра infl позволяет предотвратить развитие расходимости процесса усвоения, при этом адаптивное задание infl позволяет получить более точную оценку как прогнозируемой переменной так и неизвестного параметра.

Заключение

В докладе рассмотрены современные постановки задачи усвоения данных, основанные на ансамблевом подходе. Для повышения эффективности этих алгоритмов предлагается использовать алгоритм с трансформацией ковариационных матриц, в частности, предложенный в [1,10] ансамблевый пи-алгоритм. Рассмотрены способы управления сходимостью ансамблевых алгоритмов.

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

[1] Климова Е.Г. Стохастический ансамблевый фильтр Калмана с трансформацией ансамбля возмущений // Сибирский журнал вычислительной математики /РАН. Сибирское отделение Новосибирск. 2019. Т. 22, №1. С. 27 40.

[2] Справочник по теории автоматического управления под ред. А.А. Красовского. Москва: Наука, 1987. 711 с.

[3] Carrassi A., Bocquet М., Bcrtirio L., Everiseri G. Data assimilation in the geosciences: An overview of methods, issuers and perspectives // Wiley interdisciplinary reviews: Climate Change. 2018. V. 131, Issue 5, c535, doi: 10.1002/wcc535.

[4] Evensen, G. Data assimilation. The ensemble Kalnian filter, Spriger-Verlag, Berlin Heideberg, 2009.

[5] Evensen G. Analysis of iterative ensemble smoother for solving inverse problems. Computational Geosciences. 2018. 22. P.885-908. doi.org/10.1007/sl0596-018-9731-y.

[6] Evensen, G., P.J. van Leeuwen An ensemble Kalman smoother for nonlinear dynamics // Monthly Weather Review. 2000. V. 128. P. 1852-1867.

[7] Fisher M., Leutbecher M., Kelly G.A. On the equivalence between Kalman smoothing and weak-constraint four-dimensional variational data assimilation// Quarterly Journal of the Royal Meteorological Society. 2005. V. 131, P. 3235-3246.

[8] Houtekamer, H.L. Zhang, F. Review of the ensemble Kalman filter for atmospheric data assimilation // Monthly Weather Review. 2016. V. 144, P. 4489-4532.

[9] Jazwinski, A.H. Stochastic processes and filtering theory. Academic Press, New York, 1970.

[10] Klimova E. A suboptimal data assimilation algorithm based on the ensemble Kalman filter // Quarterly-Journal of the Royal Meteorological Society. 2012. V. 138, P. 2079-2085.

[11] Klimova E.G. Application of ensemble Kalman filter in environment data assimilation // IOP Conference Series: Earth and Environmental Science. - Vol. 211 (??) 012049 . - doi:10.1088/1755-1315/211/l/012049.

[12] Lei, J., Bickel, P., Shyder, C. Comparison of ensemble Kalman filters under non-gaussianity //Monthly-Weather Review. 2010. V. 138. P. 1293-1306.

[13] Nakamura G., Potthast R. Inverse Modeling. IOP Publishing. 2015. 509 pp. doi.org/10.1088/978-0-7503-1218-9.

Климова Екатерина Георгиевна — д.ф.-м.н., ст. науч. сотр. Института

вычислительных технологий СО РАН;

e-mail: klimova@ict.nsc.ru Дата поступления — 30 апреля 2019 г.