Оценка параметров состояния атмосферы с использованием четырехмерной динамико-стохастической модели и линейного фильтра Калмана. Часть 1. Методические основы Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

УДК 551. 501; 512.25 Ю.Б. Попов

Оценка параметров состояния атмосферы с использованием четырехмерной динамико-стохастической модели и линейного фильтра Калмана. Часть 1. Методические основы

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

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

Введение

Для решения прикладных оперативных задач в области метеорологии, геофизики, экологического мониторинга, управления воздушным движением, при чрезвычайных ситуациях и катастрофах природного и техногенного характера, а также для ряда задач оборонного значения существует потребность в текущей и прогностической информации о пространственно-временной структуре полей метеорологических величин. Большинство из перечисленных задач решаются в рамках ограниченных территорий с горизонтальными размерами 50-500 км и высотой верхней границы 10 км.

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

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

Рост требований к точности текущей и прогностической информации, а также к ее оперативности, пространственному и временному разрешению, при минимальных исходных данных, привели к необходимости разработки новых более совершенных методов. Так, в последние годы интенсивное развитие получили численные методы пространственной и временной экстраполяции полей метеорологических величин, основанные на использовании аппарата фильтрации Калмана [1, 6, 7]. Привлекательность этих методов обусловлена возможностью использования хорошо разработанных математических моделей на основе дифференциальных уравнений, в том числе нелинейных и стохастических, описывающих поведение распределенных динамических объектов в пространстве и во времени. Кроме того, эти методы позволяют синтезировать рекуррентные алгоритмы, удобные для реализации на современной микропроцессорной технике.

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

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

1. Физическая постановка задачи экстраполяции полей метеорологических величин с обоснованием моделей состояний и наблюдений

Начнем постановку задачи с предположения, что интересующие нас метеорологические величины (например, давление, температура, влажность, ортогональные составляющие скорости ветра) непрерывно распределены в некотором объеме пространства. Допустим, что этот объем ограничен снизу уровнем земли, сверху - высотой рассматриваемого слоя атмосферы h=10 км. Вид нижнего сечения определяется выбранным мезомас-штабным полигоном, максимальные линейные размеры которого не превышают (500x500) км. В пределах подобного полигона произвольно размещено 8 аэрологических станций, обеспечивающих регулярные измерения метеорологических величин во всем интервале высот до верхней границы h с заданным разрешением по высоте Ahj (/ = 1, 2, ..., Щ. Для некоторого фиксированного момента времени tf¡ данные от каждой станции можно представить в виде ^мерного профиля (вектора), каждая компонента которого является значением метеорологической величины, измеренной на соответствующем высотном уровне h = h2, hз, ..., hN.

Физическая постановка задачи пространственной экстраполяции (интерполяции) состоит в том, чтобы в текущий момент времени tk по совокупности данных от 8 разнесенных измерительных станций сформировать оценку метеорологической величины в произвольной точке заданного пространства (или узле регулярной сетки) с координатами (х0, у0,

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

Для начала определим несколько условий и введем основные обозначения. Так, для каждой аэрологической станции выберем из ^мерного профиля три значения поля одновременно измеренных на разных высотных уровнях. Первый из них соответствует высотному уровню прогнозируемой точки пространства (^ =hj, ] = 1, ...,Щ. Второй и третий - соответственно, ближайший нижний hj-l и верхний hj+l высотные уровни по отношению к первому. Также будем считать, что результирующее значения поля §0 в точке экстраполяции с координатами (х0, У0, определяется суммой регулярной ^о и флуктуа-ционной ^0 составляющих:

При этом в качестве регулярной составляющей поля в точке экстраполяции используем средневзвешенное значение для трёх ближайших станций [4], которое рассчитывается в плоскости для заданного высотного уровня ho:

3

^0 = % + ^0 .

(1)

Е ъ ■ ^

(2)

где ^ - измеренное значение поля в г-й точке; qi = 1-

весовой коэффициент;

Е р >о

13 =1

- проекция расстояния между г-й станцией и точкой экстраполяции на горизонтальную плоскость; (х, у) - прямоугольные координаты.

Соответственно, флуктуационная (случайная) составляющая поля может быть найдена из (1) следующим образом:

% = - 1о. (3)

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

ц(Дт) = ехр(-а-Дх); ц(Др) = ехр(-р-Др); ц(Дй) = ехр(-у-Дй), (4)

где ц(Дт),|а(Др),|а(Дй) - соответствующие корреляционные функции; а= 1Д , р =

= 1/ Р = V Дто /Дро

у = y.i а = у. , р = у. , у = укг> - коэффициенты, обратно пропорциональные / ДНо / ДРо / ДРо / ДНо

интервалам временной (Дро), пространственной (Дро ) и межуровневой (ДНо) корреляции

соответственно. Согласно результатам, опубликованным в [2, 8], для типового мезомас-

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

Дро =3оо км, ДНо = 45оо м (для температуры); Дто =24 ч, Дро =2оо км, ДАо =15оо м (для ортогональных составляющих скорости ветра). Очевидно, что в условиях различных ме-зомасштабных процессов эти параметры изменяются. Соответственно, неопределенность значений Дто , Дро и ДНо является ограничением в использовании (4). Для снятия подобного ограничения могут быть применены различные подходы, позволяющие выполнить оценку а, Р, у в реальном масштабе времени, с использованием текущих измерений метеовеличин. В частности, в [2] предложен подход на основе фильтра Калмана, обеспечивающий одновременную оценку искомых метеорологических величин и коэффициентов а и р. Аналогично можно выполнить оценку у. Вышеизложенное позволяет сделать предварительное допущение, что коэффициенты а, Р, у известны и постоянны на

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

Далее используем подход, предложенный автором в [3].

Используя временную корреляционную функцию (4), запишем в разностной форме для точки прогнозирования (i = о) динамико-стохастическую модель эволюции значений флуктуационной компоненты поля ^о во времени:

% (k + 1) = & (k) • (1 - а • Д) + шо (k), (5)

где At — интервал временной дискретизации; k = о, 1, ..., K — дискретное время (tk = kAt);

2

Юо^) - порождающий шум с дисперсией о , определяющий случайный характер про-

шо

цесса +1) [5].

Аналогично, используя выражение для пространственной корреляционной функции ц(Др) = exp( - в •Др), введем модель, определяющую зависимость между двумя точками

поля в горизонтальной плоскости. Запишем в разностной форме динамико-стохасти-ческую модель, связывающую значения флуктуационной компоненты поля для произвольной точки i и точки прогнозирования о:

(k + 1) = &(k+1) • (1-в •Дрю) + ш (k+1), (6)

где Др^о - проекция расстояния между i-й точкой и точкой экстраполяции на горизонтальную плоскость; Oi (k + 1) - порождающие шумы, определяющие случайный характер процессов (k +1) [5]. Физическую трактовку последнего уравнения можно представить следующим образом [2]. Если в момент времени t = Ü в некоторой точке произошло возмущение поля то в любой момент времени на расстояниях, сравнимых с Др^, отклик на это возмущение будет удовлетворять соотношению (6).

Исследования, опубликованные в [8], дают основания утверждать о наличии статистической связи между значениями метеорологических величин, измеренных на двух ближайших высотных уровнях. В этом случае из выражения ц(ДА) = exp( - у •ДА) следует

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

,йо +АН/о (к +1) = 5*,^ +1) • (1 - у • АН/о) + ® * (к +1), (7)

где АН/о - разность высот между выбранными уровнями; / = 0, 1, 2 - номер высотного уровня (0 - уровень точки прогноза, 1 - нижележащий, 2 - вышележащий уровень).

Физическая трактовка уравнения (7) аналогична (6), только в данном случае пространственным параметром будет высота.

Далее, выполняя подстановку значения Ю(к +1), из выражения (5) в (6), получим

(к +1) = ^(к) • (1 -а-АО • (1 -Р-Арю) + го* (к). (8)

Правую часть выражения (8) подставим в правую часть (7), т.е. выполним замену ^^(к +1) на (к +1). После подстановки получим окончательное выражение для дина-

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

% (к +1) = • (1 - а • АЬ) • (1 -р-Ар*о) • (1 - У • АН/о) + го* (к), (9)

где го* (К) - эквивалентные порождающие шумы, определяющие случайный характер процесса.

Порождающие шумы считаются белыми, нормальными, с нулевым средним и известной корреляционной матрицей:

М{Щк)}= о и М{П(к)-ОГ(1)}=Яа(к)-5к1, (Ю)

Т I I

где □ (к) =гоо(&),го1(к),го2(к),...,го*(к),...,го'п(к) - вектор-столбец эквивалентных шумов состояния размерностью (п+1); М{-} - оператор усреднения; 5кг - символ Кронекера; Т - оператор транспонирования;

Яо,о Яо,1 ''' Яо,п

Ко (к) =

Я1,о Я1,1 Я1,п

корреляционная матрица эквивалентных шумов состояния.

Яп,о Яп,1 Яп,п

В качестве примера приведем типичные выражения для расчета элементов матрицы Кп(к), для точки прогноза и одной станции измерения (для трех высотных уровней). При этом используем следующее выражение Я*,/ = М{го*-го/} (*, / = о, 1, 2, 3):

Во,о =<;

Я1,1 =<а -РАР1,о)2 +<;

Я2,2 =°гоо(1 -РАр1,о)2(1 -уАН1,о)2 +^(1 -уАН1,о)2 +а;;2; Rз,з = (1 -рАр10)2(1 - уАЛ2,0)2 + а;21 (1 - уА^0)2 + а2з; Яо,1 (к) = Я1,о(к) = гоо (1 -РАр1,о );

Яо,2 (к) = Я2,о(к) = а|ТОо(1 -рАр1,о) • (1 -уАН^о);

Яо,з(к) = Яз,о(к) = оТОо(1 -рАр1,о) • (1 -уАН2,о); (11)

Я1,2 (к) = Я2,1 (к) = ^гоо(1 - РАр1,о )2(1 - уАН1,о) + стго1(1 - уАА^о);

Я1,з (к) = Яз,1 (к) = а;ТОо(1 - РАр1,о)2(1 - уАН2,о) + ст^^Ц - уАН2,о);

Я2,з(к) = Яз,2 (к) = а|ТОо(1 -хАр1,о)2(1 -уАН1,о)(1 -уАН2,о) + аТО1(1 -уАН1,о)(1 -уДЙ2,о),

2 2 2 2

где ст , а , а , ст - дисперсии исходных шумов состояния для точки экстраполя-

гоо го1 го2 гоз

ции и точки измерений на трех высотных уровнях соответственно.

Расчет элементов матрицы Кд(к) для большего числа точек измерений выполняется аналогично.

Таким образом, пространственно-временная связь между значениями поля метеорологической величины для произвольной станции * (*=1, ..., 8) и точкой экстраполяции

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

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

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

(А + 1) = ^ (А + 1) + ^ (А + 1), (12)

где (А +1) - фактические измерения, полученные в точке Ь в момент времени (А+1); е, (А +1) - ошибка измерений, нормальный белый шум с нулевым средним и известной дисперсией.

Уточним, что при записи математических моделей (5), (9), (12) априори предполагалось использование трех измерений метеовеличины на разных высотных уровнях для каждой аэрологической станции. Таким образом, для измерительной сети, состоящей из 8 пространственно разнесенных аэрологических станций, число уравнений, входящих в модель состояний (5) и (9), будет равно п = (3*$ + 1). Число уравнений, входящих в модель наблюдений (12), равно в = 3*8.

На основе уравнений (5), (9), (12) выполним синтез алгоритма, обеспечивающего оптимальную оценку текущих значений метеовеличины Ю.

2. Формализация задачи в терминах фильтрации Калмана

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

X(k +1) = Ф(А) • X(k) + Щ(А), (13)

Z(k) = Н(А) • Х(А) + Е(й), (14)

где Х(А) - вектор-столбец размерностью (п*1), включающий в себя неизвестные и подлежащие оцениванию переменные состояния динамической системы (вектор состояния); А = 0, 1, 2, ..., К - дискретное текущее время с интервалом дискретизации ДЬ (¿а = АДЬ); Ф(А) - матрица перехода для дискретной системы размерностью (п*п); Щ(А) - вектор-столбец шумов состояния, размерностью (п*1); Z(k) — вектор фактических измерений размерностью (в*1); Н(А) - матрица наблюдений, определяющая функциональную связь между истинными значениями переменных состояния и измерительными каналами системы размерностью (в*п); E(k) — вектор ошибок измерений (шум измерений), размерностью (в*1).

При этом в соответствии с [5] матричное уравнение оптимального оценивания вектора состояний Х(А) имеет вид

Л Л Л

Х(А +11А +1) = Х(А +11 А) + О(А +1) • [Ъ(А +1) - Н(А +1) • Х(А +11 А)], (15)

Л

где X7(А +1| А +1) = |х1,х2,хз,...,хп| - оценка вектора состояния на момент времени (А+1);

ЛЛ

Х(А +1| А) = Ф(А) • X(А | А) - матричное уравнение для расчета вектора предсказания оценок на момент времени (А + 1) по данным на шаге А; О (А +1) - матрица весовых коэффициентов размерностью (п*в).

Расчет весовых коэффициентов осуществляется по рекуррентным матричным уравнениям вида

О(А +1) = Р(А +11 А) • Нт (А +1) • [Н(А +1) • Р(А +11 А) • Нт (А +1) + R Е (А +1)]-1, (16) Р(А +11 А) = Ф(А) • Р(А | А) • Фт (А) + RQ (А +1), (17)

Р(А | А) = [I - О(А) • Н(А)] • Р(А +11 А), (18)

где Р(А +11 А) - апостериорная матрица ковариаций ошибок предсказания размерностью (п*п); Р(А +1| А +1) - априорная матрица ковариаций ошибок оценивания размерностью

(п*п); R_E (к +1) - диагональная матрица ковариаций шумов наблюдения размерностью

(вхв); RQ(к +1) - диагональная матрица ковариаций шумов состояния размерностью

(п*п); I - единичная матрица размерностью (п*п).

Сопоставим матричные уравнения (15)-(18) с моделью (5), (9), (12) и формализуем необходимые векторы и матрицы, входящие в уравнения фильтрации (15)-(18).

Определим вектор состояния в виде совокупности подлежащих оцениванию метеовеличин и ^:

Хт (к) = \х1(к),х2(к),хг(к),...,хп (к) =|^Ь(k),^1(k),й2(k),...,^зs (к), (19)

где Х1(к), Х2(к), Хз(к), ..., хп(к) - элементы вектора состояний, соответствующие искомым метеовеличинам.

Определим матрицу перехода для уравнений состояния в виде

(1 -а-Д^

Ф(к) =

(1 -а•Дt) - (1 -Р-Дрю) (1 - а • Д^ • (1 -Р-Дрю) • (1 -уД^2о) (1 - а • Д^ • (1 -р-Др2о) • (1 -У-Д^ю)

(1 - а • Д^ • (1 -Р^о) • (1 -У-Д^2о) 0

0 • • 0

0 • • 0

0 • • 0

0 • • 0

0 • • 0

0 • • 0

(20)

Теперь запишем матрицу наблюдений Н(к):

0 1 о ... о

Н(к) =

оо1 ооо

ооо

(21)

Таким образом, матрицы (2о) и (21) определяют структуру линейного фильтра Кал-мана (ФК) [5], обеспечивающего оценку всех элементов вектора состояний (19).

Для инициации алгоритма фильтрации (15)-(18) в момент к = о необходимо задать

л

следующие начальные условия: Х(о) = М"{Х(0)} - вектор начальных оценок состояния; Р(о|о) = М{[Х(0) - М{Х(0)}][Х(0) - М{Х(0)}]Т} - начальная матрица ковариаций ошибок оценивания, а также значения элементов ковариационных матриц шумов RE(0) и Rn(0).

л

На практике значения Х(0) и Р(0 | 0) могут быть заданы исходя из минимального объема сведений о реальных свойствах системы, а в случае полного отсутствия полезной инфор-

л

мации задаются Х(0) = 0 , а Р(0|0) = I.

В процессе работы ФК одновременно будут оцениваться неизвестные значения Ю и измеряемые с некоторой погрешностью значения ^. Итоговая оценка метеовеличины в точке экстраполяции находится по формуле

^о(к) = Х1(к) -!о(к), (22)

где Х1(к) - текущая оценка центрированной величины Ю; ^о(к) - регулярная составляющая поля в точке экстраполяции. 3. Заключение

В заключение следует отметить, что модель (5), (9), (12) позволяет наращивать число измерительных станций, а также число вовлеченных высотных уровней. Перспективы развития метода следует рассматривать в направлении предварительной или одновременной оценки условно известных параметров а, Р, у за счет расширения вектора состояния [2].

Об эффективности предложенного подхода, точности оценивания и прогноза мезоме-теорологических полей можно будет судить по данным численного эксперимента. Результаты исследований с привлечением реальных многолетних аэрологических измерений, полученных на типичном Европейском мезометеорологическом полигоне, будут являться предметом второй части этой статьи.

Литература.

1. Климова Е.Г. Методика усвоения данных метеонаблюдений на основе обобщенного субоптимального фильтра Калмана // Метеорология и гидрология. - 1997. - №11. - C. 55-65.

2. Комаров В.С. Динамико-стохастические методы и их применение в прикладной метеорологии / В.С. Комаров, Ю.Б. Попов Ю.Б., С.С. Суворов, В.А. Кураков. - Томск: Изд-во ИОА СО РАН, 2004. - 236 с.

3. Комаров В.С. Оценивание и прогнозирование параметров состояния атмосферы с помощью алгоритма фильтра Калмана. Ч. 1. Методические основы / В.С. Комаров, Ю.Б. Попов // Оптика атмосферы и океана. - 2001. - Т. 14, № 4. - С. 255-260.

4. Пространственная экстраполяция метеорологических полей в области мезомасшта-ба на основе четырехмерной смешанной динамико-стохастической модели и аппарата калмановской фильтрации / В.С. Комаров, А.В. Лавриненко, Н.Я. Ломакина, Ю.Б. Попов, А.И. Попова, С.Н. Ильин // Оптика атмосферы и океана. - 2004. - Т. 17, № 8. - С. 651-656.

5. Сейдж Э.П. Теория оценивания и ее применение в связи и управлении / Э.П. Сейдж, Дж.Л. Мэлса. - М.: Связь, 1976. - 496 с.

6. Dee D.P. Simplification of the Kalman filter for meteorological data assimilation // Q.J. Roy. Met. Soc. - 1991. - Vol. 117. - P. 365-384.

7. Ide K. Unified Notation for Data Assimilation: Operational, Sequential and Variation / K. Ide, Ph. Courter, M. Ghil, A. Lorenc // Jour. Of Meteorological Society of Japan. -1997. - Vol. 75, № 1B. - Р. 1120-1126.

8. Комаров В.С. Методика сверхкраткосрочного прогноза параметров состояния атмосферы на основе алгоритма калмановской фильтрации и двумерной динамико-стохастичес-кой модели / В.С. Комаров, А.В. Лавриненко, Ю.Б. Попов // Оптика атмосферы и океана. - 2005. - Т. 18, № 4. - С. 344-348.

Попов Юрий Борисович

Канд. техн. наук, доцент каф. радиоэлектроники, Сургутского гос. университета

Тел.: 8-922-405-05-47

Эл. почта: popovyub@mail.ru

Popov Yu.B.

The estimation of atmospheric state parameters using a four-dimensional dynamic-stochastic model and a linear Kalman Filter. Part 1. Methodical basis

In the article there are discussed the questions of algorithm synthesis for parameters of spatial atmospheric state forecast within the local territories (e.g. of temperature and wind). The algorithm is based on Kalman filtering techniques and a new four-dimensional low-parametric dynamic-stochastic model. The results are of interest for solving the applied problems in such fields as meteorology, geophysics, ecological monitoring, air traffic control, at emergency situations and any accidents, and also are used for a number of problems of defense value.

Keywords: Kalman filter, spatial interpolation, data assimilation, numerical modeling, mesoscale.