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

Ранее в ИКИ РАН была разработана технология оценки запаса стволовой древесины лесов России за период с начала этого столетия. Получаемая по спутниковым данным динамика запасов стволовой древесины важна для характеристики изменения биомассы и углерода в лесах страны. Но анализ этих данных показывает, что погрешность оценки запаса стволовой древесины составляет порядка 10 % от значения запаса. Например, для леса с запасом 200–300 м3/га погрешность может составлять 20–30 м3/га. Но типичный ежегодный прирост запаса значительно ниже этой величины. Например, в таблицах хода роста (Швиденко А.З., Щепащенко Д.Г., Нильссон С., Булуй Ю.И. Таблицы и модели хода роста и продуктивности насаждений основных лесообразующих пород Северной Евразии (нормативно-справочные материалы). 2-е изд., дополн. М., 2008) даже для быстрорастущих пород в наиболее благоприятных условиях прирост не превышает 20 м3/га в год, а для типичных для России средневозрастных и более старых лиственничных насаждений III и IV бонитета их прирост составляет 1–3 м3/га в год и менее. Но при этом предполагаемые изменения за весь период спутниковых наблюдений с 2002 по 2023 гг. при предположении устойчивого роста превышает погрешность, что даёт возможность оценки долгосрочной динамики запаса стволовой древесины и, соответственно, потоков углерода. Для оценки этой динамики необходимо выполнить фильтрацию данных, отделив фактическое изменение запаса от случайного шума. На рисунке 1 показан пример временного ряда запаса стволовой древесины в одной точке, с относительно небольшим уровнем общей погрешности, который, тем не менее, превышает экспертную оценку фактической динамики запаса.

   

Рисунок 1 — Пример временного ряда запаса стволовой древесины с погрешностью на уровне 10 % (270±30 м3/га), малой в абсолютном выражении (слева), но большой с точки зрения предположительной фактической динамики запаса (справа)

         
а)                                                                               б)
         
в)                                                                            г)
Рисунок 2 — Пример динамики запаса стволовой древесины по данным ДЗЗ и отклонения результата его сглаживания от предположительной фактической динамики. Для рисунка г не приведена фактическая динамика, так как она совпадает с графиком STL

На данный момент для решения задачи уточнения динамики запаса стволовой древесины в ИКИ РАН используется алгоритм сглаживания STL (Cleveland R.B., Cleveland W.S., McRae J.E., Terpenning I. 1990. STL: A Seasonal-Trend Decomposition Procedure Based on LOESS // J. Official Statistics. V. 6. P. 3–73). Он раскладывает временной ряд на трендовую, сезонную и случайную составляющие, используя для построения тренда локальную взвешенную регрессию (алгоритм LOWESS). Данный метод показывает достаточно хорошие результаты для значительной части лесов страны, но может приводить к неверной оценке динамики запаса в ряде случаев. Это демонстрирует рисунок 2, на котором показаны фактические оценки запаса по данным ДЗЗ, результат сглаживания алгоритмом STL, а также экспертная оценка вероятной фактической динамики запаса стволовой древесины для каждого из четыре рассмотренных случаев. Например, на рисунке 2а показан пример результата работы алгоритма STL, демонстрирующий сложную динамику с ростом и падением запаса, что противоречит экспертной оценке динамики стабильного медленного роста в старовозрастных лесах. На рисунке 2б показан пример наличия больших отклонений значений запаса в первые несколько лет, что приводит к ошибочной положительной динамике после применения алгоритма STL. При этом на рисунках 2в, г сглаживание STL адекватно оценивает динамику запаса при тех же настройках метода. Стоит отметить, что рисунок 2г демонстрирует пример изменения скорости роста в течение 20 лет наблюдения, что не позволяет просто аппроксимировать ряд линейной функцией. Большое разнообразие возможной динамики запаса стволовой древесины приводит к необходимости разработки адаптивного алгоритма, способного адекватно описывать динамику запаса для любых лесов на территории страны. В отчётный период был опробован подход к оценке фактической динамика запаса стволовой древесины на основе фильтра Калмана (Welch G., Bishop G. An Introduction to the Kalman Filter. 1995). Фильтр Калмана предназначен для оценки фактического состояния системы при наличии её модельного описания и наблюдений её состояния, обладающих известными погрешностями.

Базовый фильтр Калмана предполагает линейную модель изменения системы

где x_k– состояние системы в момент k, F – линейное преобразование, описывающее эволюцию системы, ω_k — случайные шум, распределённый по нормальному закону ω_k~N(0,Q_k), где Q_k– матрица ковариации.

Аналогично формула, описывающая наблюдения состояния системы с погрешностями, имеет вид

где z_k – наблюдения состояния системы, H – линейное преобразование между состоянием системы и фактическими наблюдениями, и v_k~N(0,R_k) – случайный шум наблюдений, распределенный по нормальному закону.

При наличии такого описания системы фильтр Калмана позволяет найти оптимальную оценку состояния системы с учётом погрешности модели и наблюдений. На первой этапе фильтр Калмана оценивает априорные значения состояния системы только по модели, без учёта наблюдений, задаваемые формулами

где – априорная оценка состояния в момент времени k, P_{k | k-1} – априорная матрица ковариации ошибки в момент времени k.

Затем происходит оценка оптимального состояния системы с учётом наличия наблюдений, для чего оценивается отклонение фактических наблюдений от ожидаемой (априорной) модельной оценки

После чего вычисляются необходимые для определения оптимального состояния системы матрица Калмана в момент времени k – S_k, и коэффициент усиления Калмана в момент времени k – K_k:

На основе которых оценивается оптимальное (апостериорное) состояния системы () и обновленная оценка матрицы ошибки (P_{k | k}):

Ввиду того, что задача может быть сформулирована как сглаживание временных рядов запаса, была использована сглаживающая версия фильтра Калмана (Kalman smoother), отличающаяся от стандартной версии наличием обратного прохода по временному ряду. Для выполнения сглаживания использовался метод Рауха – Тунга – Штрибеля (Rauch H.E., Tung F., Striebel C.T. Maximum likelihood estimates of linear dynamic systems // AIAA J. 1965. V. 3.No. 8. P. 1445–1450), с обратным проходом по временному ряду от k = n–1 до 1, выполняемого на основе следующих формул:

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

Для наблюдений используется следующая формула:

с использованием фактических спутниковых наблюдений GSV_k^RS, и рассчитанного по ним наклона линейной регрессии LINEAR.
Для оценки матриц погрешности модели Q и наблюдений R используются следующие формулы:

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

где GSV_ERROR – погрешность данных ДЗЗ, полученная ранее в рамках метода их построения, а LINEAR_ERROR – погрешность коэффициента наклона линейной регрессии.

Приведённые выше формулы были реализованы на языке Python и применены к временным рядам оценок запаса стволовой древесины по данным ДЗЗ. Результат их применения приведён на рисунке 3. На каждом из графиков видно, что фильтр Калмана с описанными выше параметрами может адекватно описывать динамику запаса стволовой древесины, компенсируя недостатки используемого ранее метода сглаживания на основе алгоритма STL. На рисунке 3а фильтр Калмана показывает устойчивую линейную динамику, хорошо совпадающую с экспертной оценокй фактической динамики, без периодов роста и падения запаса, демонстрируемых методом STL. На рисунке 3б фильтра Калмана также значительно лучше воспроизводит экспертную оценку динамики, игнорируя сбойные выбросы в 2004-2006 годах. Также фильтра Калмана показывает адекватную динамику запаса стволовой древесины на рисунке 3в, и способен воспроизвести изменение темпов прироста запаса на рисунке 3г.
По результатам проведённого анализа можно заключить, что фильтр Калмана может применяться для сглаживания временного ряда запаса стволовой древесины и оценки его динамики. В дальнейшем предложенный метод на основе фильтра Калмана будет применён для всех лесов России. Для этого потребуется его модификация, учитывающая возможность внезапных изменений значений запаса стволовой древесины из-за нарушений, таких как пожары и рубки. Будет рассмотрена возможность интеграции в фильтр Калмана дополнительных спутниковых наблюдений, детектирующих факт таких нарушений.

       
а)                                                                                  б)
         
в)                                                                              г)
Рисунок 3 — Пример применение фильтра Калмана к динамике запаса стволовой древесины лесов по данным ДЗЗ в сравнении с полученными ранее результатами на основе метода STL