CC BY
63
2006
Известия ТИНРО
Том 145
УДК 573.22.087.1.001.57
А.А.Михеев (СахНИРО, г. Южно-Сахалинск)
ПРИМЕНЕНИЕ ОБОБЩЕННОГО ФИЛЬТРА КАЛМАНА К ПРОГНОЗИРОВАНИЮ ВОЗВРАТОВ ГОРБУШИ*
Представлен новый алгоритм оптимизации параметров для обобщенного фильтра Калмана 1-го порядка, базирующийся на аналитическом вычислении градиентов логарифмического правдоподобия. Проведенные исследования позволили обосновать сходимость указанного фильтра для модели учета одновозрастной популяции, описываемой функцией Риккера. Полученные результаты анализа были подтверждены при моделировании динамики нерестовых возвратов горбуши юго-восточного Сахалина за период 1970-2002 гг. С помощью построенного обобщенного фильтра Калмана реконструировали динамику подходов названного промыслового стада и дали ее прогноз на 2004 г.
Mikheyev А.А. Application of the Extended Kalman Filter for forecasting of pink salmon returns // Izv. TINRO. — 2006. — Vol. 145. — P. 146-167.
A new algorithm of parameter optimization for the Extended Kalman Filter (EKF) of the first order is presented. It bases on analytical calculation of the log likelihood gradient. Convergence of the filter is proved for accounting of the same-aged population described by the Ricker function. The obtained results are confirmed by the modeling of spawning returns for the southeast Sakhalin population of pink salmon in 1970-2002. The return in 2004 was forecasted.
Современное управление морскими биоресурсами построено на определении допустимого изъятия на некоторую перспективу. Понятно, что обоснованно установить необходимый уровень изъятия без анализа динамики самих ресурсов невозможно. Однако размеры промысловых популяций, как правило, недоступны для прямого наблюдения. Исходная информация для анализа обычно состоит из косвенных характеристик: учетных и промысловых индексов численности, годовых уловов, выборочного размерного состава и навесок. Если расчеты не отражают связей запаса с измеренными характеристиками, то в полученной оценке возникает методикообусловленная неопределенность (Francis, Shotton, 1997; Chen, Wilson, 2002). Другим неустранимым источником неопределенности являются погрешности самого измерения. Они вызваны искажениями и неполнотой промысловой статистики, а также неизбежными случайными ошибками прямого учета, связанными с определением размеров обловленного пространства, влиянием среды на уловистость орудия лова и т.д. Кроме того, на состояние запаса в краткосрочном масштабе воздействует комбинация факторов внешней среды, априори неизвестная и практически непредсказуемая. В результате полученные измерения несут помимо собственных ошибок еще и следы шума, генерируемого при-
* Публикуется в дискуссионном порядке. Работы в развитие этой темы будут публиковаться вне очереди. — Ред.
родными процессами. Резюмируя, можно сказать, что системная динамика, т.е. динамика величины запаса, и связанная с ней последовательность измерений являются двумя разными, но сопряженными стохастическими процессами. Отказ принимать во внимание данное обстоятельство, возможно, одна из главных причин неудачного регулирования промысла, что на самом деле довольно распространенное в мировой практике явление (Schnute, Richards, 2001). Таким образом, в сфере принятия решений по управлению морскими биологическими ресурсами, в той ее части, которая касается определения допустимых уловов, существует проблема методического характера, состоящая в игнорировании структурных свойств и способа получения исходной информации. Указанная проблема тем более актуальна, что значительная часть ресурсов находится сегодня в подавленном состоянии, а объем и качество поступающей информации продолжают снижаться ^орельский, Макоедов, 2004).
Обозначенная выше проблема является предметом теории фильтрации и сводится к задаче удаления шума и выделения полезного сигнала из входного информационного потока (Вероятность и математическая статистика, 1999). Указанная теория восходит к методу наименьших квадратов (MHK) (Sorenson, 1970). Стимулированная потребностями в надежных системах связи, слежения и ориентации, она получила существенное развитие в работах А.Н^олмогорова (1941) и Винера (Wiener, 1949). Однако применение числовых фильтров в прикладных исследованиях и технике стало реальностью только после появления работ ^лмана и Бьюси (Kalman, 1960; Kalman, Bucy, 1961). Фильтр ^лмана (ФЮ был построен как оптимальный рекурсивный алгоритм обработки данных и впервые реализован в системе ориентации космического корабля Аполлон-XI и его лунного модуля. В настоящее время ФK встречается в самых разнообразных областях, это: космическая и наземная навигация, искусственный интеллект и обучающие нейронные сети, автоматическое регулирование и управление на основе промышленных АСУ, роботы, связь и цифровая обработка сигналов, мониторинг и прогноз природной среды, эконометрика и анализ фондовых рынков и ряд других.
Одним из первых на необходимость стохастических моделей с разделением шумов системы и наблюдения в рыбохозяйственных исследованиях указал Шнюте (Schnute, 1991), он представил ФK на примере простой продукционной модели. С того времени данный фильтр использовали во всех базовых моделях рыболовства (Sullivan, 1992; Pella, 1993; Gudmundsson, 1994, 1998, 2004; Schnute, 1994; Kimura et al., 1996; Reed, Simons, 1996). ^уг задач, решаемых с его помощью, весьма широк. Это и определение биологических ориентиров на основе моделей "запас—пополнение", и анализ влияния внешней среды на продуктивность ряда запасов тихоокеанских лососей (Peterman et al., 2000, 2003; Schnute, Kronlund, 2002). В последних работах динамику популяций описывали функцией Риккера, которая, как известно, является базовой моделью для отражения связи обилия родителей и потомков, особенно у лососей (Hilborn, Walters, 1992). Названная функция также была использована в ФK и автором — при прогнозировании нерестовых подходов горбуши в два района о. Сахалин Михеев, 2003, 2004а, б).
Обработка данных на основе ФK в рыбопромысловых моделях неразрывно связана с параметрической идентификацией, т.е. с определением констант, при которых достигается лучшее соответствие между модельной динамикой и наблюдениями (Schnute, 1991). Если динамика системы описывается нелинейной функцией, то приходится применять обобщенный ФK (ОФЮ (Справочник по прикладной статистике, 1990). Отсюда возникает задача поиска глобального максимума у нелинейной целевой функции, в качестве которой обычно берется правдоподобие прогнозной оценки. Поиск, как правило, осуществляется итерационными градиентными методами (Банди, 1988). ^гда обновление информации представляет собой стационарный гауссовский процесс, применяются специальные градиентные методы, относящиеся к классу нелинейных MHK (Демиденко, 1981). Такая
ситуация типична для традиционных областей применения ФК (Справочник по прикладной статистике, 1990).
Если же ковариация обновляющего процесса зависит от параметров модели, то требуется привлечение градиентных методов общего типа. Для изменяющейся во времени ковариации приходится определять ее значение и градиент при каждом новом поступлении информации, что без соответствующих рекуррентных формул становится чрезвычайно сложной вычислительной задачей (Справочник по прикладной статистике, 1990). В доступной автору литературе, значительная часть которой перечислена в этой работе, указанная задача не рассматривалась. В связи с этим цели данной работы — получить требуемые рекуррентные формулы и на их основе построить алгоритм для параметрического оценивания ОФК, а также продемонстрировать работу этого алгоритма при прогнозировании динамики одновозрастной популяции, подчиняющейся функции Риккера. В качестве объекта приложения выбрана тихоокеанская горбуша Oncorhynchus gorbuscha. Сделанный выбор во многом обусловлен чрезвычайной актуальностью прогнозирования этого промыслового вида.
1. Представление дискретного ОФК
Пытаясь определить состояния одного случайного процесса по значениям другого, связанного с ним, случайного процесса, приходим к стандартной постановке задачи фильтрации для модели скрытых состояний марковского процесса (рис. 1).
yt-2 yt-1 У»
w w
p(y»lx»)
Измерения
Ненаблюдаемые состояния
P(x/lx/-i)
Рис. 1. Направленный ацикличный граф для скрытых состояний марковского процесса (по: Van der Merwe, 2004; с небольшими изменениями)
Fig. 1. Oriented acyclic graph for latent states of the Markov's process (by: Van der Merwe, 2004 with some changes)
Скрытые состояния можно рассматривать для любого момента времени t как векторную переменную xt = (xt .., x ), представляющую, например, ненаблюдаемые численности m возрастных классов в популяции. Тогда временная
последовательность n точек xt е Rm, t = 1,...,n образует в пространстве (фазовом) состояний Rm дискретную динамическую систему. Пусть в каждый дискретный момент t поступает информация в виде измерения y t е Rk некоторых k характеристик данной системы (рис. 1). Тогда модель скрытых состояний марковского процесса можно представить как динамическую модель пространства состояний (Van der Merwe, 2004). Как известно, такая модель состоит из двух уравнений: системы (процесса) и наблюдения (Справочник по прикладной статистике, 1990):
x t = f (Vi,u t,v г; 6);
y t = h(xt,w t ;6),
(1.1) (1.2)
где V — шум системы (процесса), V{ е Rт; w — шум (погрешность) измерений, w{ е Rк; — известная переменная внешнего управляющего сигнала, и{ е R'; 0 — параметры, 0 е Rр; I и h — дифференцируемые функции со значениями соответственно из и Ы4.
С помощью модели (1.1), (1.2) для марковского процесса скрытых состояний могут быть определены функция априорной плотности переходной вероятности* р(х^|х ) и правдоподобие р(у^х^) следующим образом:
iх¿-1 ) = 18(Х - f (,и,,v,;0 , ^,, ру ¿1х ы8 (у *- ь(х'w *;0 *)dw *,
где 5(*) — дельта-функция Дирака; р(^) и p(wt) — известные плотности шума. Оптимальные оценки переменных X для модели (1.1), (1.2) определяются
t-х t
как условное среднее марковского процесса:
х t = e((| yiKt )=j x tp( |yl..t )xt,
где yj = (yp...,yt) — составной вектор; e(-) — оператор математического ожидания.
Отфильтрованная апостериорная плотность p(xt|yt t) может быть найдена в рамках общей теории рекурсивного байесовского вероятностного вывода (Welch, Bishop, 2001; Van der Merwe, 2004) согласно известному правилу Байеса:
( I ) p(ytlxt)pxtк.,-1)
)= f ( I w I-. (13)
j p(y t Iх t )(xt |yiK.t-1)x t
Априорную плотность p(xt|yt t-1) задает модель вероятностного процесса:
pxt |yl„t-i )=j pxt iх t-i)(x t-i |yl„t-i)x t-i, (1.4)
которая и замыкает рассматриваемый рекурсивный алгоритм.
Калман (Kalman, 1960) представил оценку Xt как взвешенную сумму текущего измерения y и априорной оценки (прогноза) x+:
x t = K' x + + K t y t. (1.5)
Тем самым он получил линейную оценку со свойствами рекурсивности. Такой выбор позволил рассматривать апостериорную плотность p(xt|yt t) в формуле Байеса как линейную комбинацию априорной плотности и правдоподобия. В результате интегральную модель вероятностного процесса, описываемую формулой (1.4), удалось свести к линейным уравнениям для средних и ковариаций соответствующих плотностей. Были определены точные аналитические значения весовых матриц K' и Kt, при которых полученная оценка становилась эффективной, т.е. была несмещенной и имела минимальную дисперсию в своем классе (Крамер, 1975).
В случае линейной модели с белым шумом названный автор получил элегантное аналитическое замкнутое решение, которое и известно как ФК, или теорема Калмана (Волков и др., 2003). Пусть уравнения системы и наблюдения имеют вид:
* Функции плотности вероятностей далее будем именовать просто плотностями.
xt = ft-ix t-1 + G t u t + V t;
(1.6)
: HX + W,
(1.7)
где Ft(0) — переходная матрица состояний размерности тхт; Gt(0) — матрица управлений размерности тх1; Н^(0) — матрица измерений размерности kхm; V и w — переменные белого шума с нулевым вектором средних и ковариационными матрицами Q и Ы размерности соответственно тхт и kхm:
е(у ) = 0; е^ )= 0; е(у (, wTг )= 0; е( , уТ )= Qдh; е(, wT )= R8tz,
где 8Н — дискретная дельта-функция Дирака (символ Кронекера), Т — символ транспонирования.
Тогда искомое решение имеет вид:
x ( =[1 - К (Н (] ++ К (у (,
где I — единичная матрица. Вес К известен в сегодняшней литературе как коэффициент усиления Калмана (Справочник по прикладной статистике, 1990). Последовательность
У» = У» - н (х;
можно рассматривать как обновляющий процесс, корректирующий априорную
оценку, и выразить через нее полученное решение: х» = х ++ к»~.
Полный алгоритм ФК для линейной модели (1.6), (1.7) показан на рис. 2.
Начальные и априорные данные
X 0, P0, Ft, G t, H t, u t, Q, R
Предиктор
1. Прогноз переменной состояния:
x t+i = Ftx t + G t+iu i+1.
2. Прогноз ковариации ошибки:
pt+i = ft p FT + q.
/V
t=t+i С
Корректор
1. Вычисление коэффициента усиления:
k t = p+hf [h t p+hf +
r]-1
2. Обновление оценки:
x t = x ++ K t ~ t •
3. Обновление ковариации ошибки:
Pt =[I - IK t H t ]P+
Текущие оценки X t
Текущее измерение yt
Рис. 2. Алгоритм типа "предиктор-корректор" для дискретного линейного ФК (по: Welch, Bishop, 200l; с небольшими изменениями)
Fig. 2. Algorithm in the form of "predictor-corrector" for the discrete linear FK (by: Welch, Bishop, 2001 with some changes)
В общей модели (1.1), (1.2) среднее и ковариация плотностей нелинейно зависят от оценки X t, и метод Калмана прямо к ним не применим. Однако функции f и h в этой модели можно представить в окрестности оценки с помощью разложения Тейлора. Если отбросить члены (п+1)-того порядка, то к аппроксимированной системе можно применить технику рекурсивного вывода. Этот подход получил название ФК n-того порядка. Для нелинейной системы с аддитивным управляющим сигналом и таким же стационарным шумом:
Xt = f (x t-i; 0)+ u t + V t; (Ш
150
y t = h(xt; 0)+ w t, (1.9)
можно построить ФК 1-го порядка со следующим алгоритмом "предиктор-
корректор :
x++i = ft + u
t+1
(1.10) (1.11)
(1.12)
(1.13)
(1.14)
Р+ = [ л ] [V xf ,]т + q; к , = Р,+ [ xh , ] [[V ,11 , ]pt+[v xh , ] + RI1;
X, = X ++ К ,[у , - h ,];
Р =У - К , [V X h , ]]+
где ft = f(х,;0); h, = ^х+;0); Vх — оператор градиента в пространстве состояний.
Для реализации представленного алгоритма остается задать параметры 0, управление и^, ковариации шума Q и Ы, а также начальные условия:
x0 = е[х0]; р0 = е[(х0 -б[х0]Хххо -е[х0])].
Именно этот алгоритм и будет использован далее как основа для моделирования.
2. Модель учета численности нерестовой горбуши
Одновозрастной состав стада горбуши позволяет построить простую модель для переменной с одним состоянием. Уравнение системы будет скалярным нелинейным:
IF (nt; a, b)+ ut+i + st, если > nm
t = 0,1,...,n. (2.1)
Учет подошедшей горбуши можно выразить простым уравнением наблюдения:
У = Щ + 8,. (2.2)
В формулах (2.1), (2.2) п — численность подхода; у — учтенная численность; ^(п,;а,Ь) — функция Риккера, описывающая связь родительского и дочер-
t
него поколений:
F(nt;a,b) = nt exp(a-bnt),
(2.3)
где a и b — положительные параметры, такие, что exp(a) отражает репродуктивный потенциал (для эксплуатируемого стада с учетом промысловой гибели), а (1/b) — емкость среды, обусловленная ограничением жизненных ресурсов (Ска-лецкая и др., 1979).
Параметром nmin в уравнении (2.1) характеризуется нижний порог численности популяции, обусловленный депенсаторным эффектом Олли (Chen et al., 2002). В уравнениях (2.1) и (2.2) величины е и Sf — переменные гауссовского белого шума с нулевым средним и дисперсиями op2 и <rm2, характеризующими соответственно интенсивность шума системы и погрешность измерений.
Аддитивный управляющий сигнал u в уравнении (2.1) полагается известным и может быть вычислен по формуле:
ut =Е Akcos
2* t
+ Vk
(2.4)
где Ak, Tk, и — параметры, соответствующие амплитуде, периоду и фазе k-той гармоники, со значениями из табл. 1. Переменная u трактуется как результат влияния внешней среды на системную динамику (Михеев, 2004а).
Таблица 1
Параметры квазигармонического тренда остатков модели Риккера для горбуши юго-востока о. Сахалин за 1970-2003 гг.
Table 1
Parameters of the Ricker's model quasi-harmonic trend residuals for the even-year southeastern Sakhalin pink salmon returns in 1970-2003
Компонента Параметр
тренда, k T., лет k' A,, млн экз. k' ф4> рад-
1 2,000 4,2-1014 1,571
2 3,362 4,948 2,466
3 5,136 1,751 2,465
4 12,098 3,906 1,618
5 34,118 6,599 1,800
Формулы (2.1 )—(2.4) вместе образуют модель учета численности половозрелой горбуши. В уравнениях (2.1)-(2.3) дискретность по времени составляет два года, а в уравнении (2.4) — один год.
Алгоритм "предиктор-корректор" ОФК 1-го порядка, соответствующий формулам (1.10)—(1.14), для модели учета горбуши примет следующий вид:
+ \Ft + ut+,, если > n min n++1 = <| t t+1' ; (2.5)
Kin. ИНаЧе
= F't2 Pt +a2p; (2.6)
P+
K = ^P4r; (2.7)
' p++* m
n = n++ Kt [ - n+]; (2.8)
Pt =(1 - Kt )P+, (2.9)
где Ft = F(nt; a,b) и Ft' = dF(; a, b)/d flt — принятые сокращения.
Чтобы реализовать представленный алгоритм, нужно задать начальные значения: среднее П0 = e[n0] и дисперсию P0 = e[(n0 — n0 )2 ]= d[n0 ]. Для рыбопромысловых моделей эти данные, как правило, неизвестны. Однако их можно определить, предположив, например, стационарность процесса в допромысловый период: e[n1 ]= e[n0] и d[n1 ] = d[n0] (Schnute, 1991; Reed, Simons, 1996). Линеаризуя
функцию F(n0;a,b) в окрестности оценки n0, легко получить уравнения для нахождения П0 и P0:
e[n1 ]= e[[ (я0; a, b)+ u1 + f1 ]= e[[ (n0; a, b )]+ u1 = F0 + u1;
d [n ] = d[ (щ; a, b)+ щ + ^ ]= fo'pd[ ]+ a P = f'p p +
152
Применение условия стационарности дает решение* no = Fo + U1>
P =^
P 1 - F'2
Другой подход состоит в произвольном выборе /?(+ при неограниченной дисперсии, т.е. в задании диффузной априорной плотности оценки (Pella, 1993). Если перегруппировать уравнения (2.7)-(2.9) для t = 0 следующим образом:
P+ +
"(y» - < );
P = P+ 1 o 1 0
(Po+))
p+
Л
+ ff 2
то при P0+ ^го они дают: Я0 = y0 и P0 = ^. Еще один способ построен на предположении, что начальное состояние измерено точно, т.е. Я0 = y0 и P0 = 0 (Schnute, Kronlund, 2002).
Заметим, что произвольный выбор начальных условий для ФК оправдан в тех случаях, когда дисперсия оценки стабилизируется: фильтр при достаточной длине ряда успевает настроиться и дает оптимальные оценки. По данному поводу имеется мнение, что логарифмическое правдоподобие слабо чувствительно к умеренной неопределенности в величине P0 (Schnute, 1994).
Однако более предпочтительным видится иной подход к проблеме выбора начальных данных: системную переменную задать, например как Я0 = y0, а дисперсию P0 считать одной из констант модели, требующих определения. При таком варианте рассматриваемая модель содержит пять неизвестных параметров:
6 ={я,b, po, ffp , ат}.
(2.10)
О соответствии модели данным у0 п = (у0,уп^ при фиксированных параметрах 0 можно судить по "концентрированному" логарифмическому правдоподобию (Ре1егтап е! а1., 2000):
L( y0K„ le) = ln \
П D (e)
t=0
exp
_ i ^ ¿(e)'
2 ¿0 Dt (e)
=_ 21Ь(eDUH»
где р{ — гауссовский обновляющий процесс (Справочник по прикладной статистике, 1990); Dt — его дисперсия.
Искомые параметры являются решением оптимизационной задачи, известной как метод максимального правдоподобия (ММП):
L(Уо...п Ö*) = maxL(Уо...„|б).
(2.12)
+
* У функции Риккера есть два отрезка, на которых > 1 и где, следовательно, значение Р0 не определено. Таким образом, этот способ выбора Р0 для функции Риккера может оказаться неприемлемым.
3. Поиск минимума целевой функции
При практическом решении задачи, представленной формулой (2.12), удобнее использовать целевую функцию, дающую тот же результат, в слегка измененном виде:
I (Уо...п |0) = -2Ь(Уо„м |в). (3.1)
Тогда задача (2.12) эквивалентна задаче поиска минимума функции (3.1):
I (Уо. К) = пнп I (Уо. |0). (3.2)
' в
Если процесс р( стационарный или быстро к нему сходится, то ММП и МНК становятся эквивалентными. Это справедливо для всех линейных систем с постоянными параметрами.
Все градиентные методы осуществляют поиск как рекуррентный процесс, представимый общей формулой (Математический энциклопедический словарь, 1988):
0 ,.+1 = 0+ А.С г.q., (3.3)
где 0. — ¿-тая оценка вектора 0; [ — номер итерации, г = 0,1,... ; qi — вектор, задающий направление поиска; G. — положительно определенная матрица, корректирующая направление поиска; — положительная константа, корректирующая шаг поиска с тем, чтобы достичь минимума на выбранном направлении. Матрица Gi и вектор q. являются функциями 0, вычисленными при значении 0 = 0..
В качестве критерия остановки вычислений обычно используют условия:
|0 .+1- 0.1 ^ |0 ,-к + *, (3.4)
где g. — градиент минимизируемой функции; £ £2 — пороговые константы, рекомендуются порядка 10-5-2-10-4 (Демиденко, 1981; Брандт, 2003); к — константа, вводимая на случай, если 0. = 0, к может быть меньше £ (Брандт, 2003). Второе условие в формуле (3.4) есть приближение необходимого условия экстремума: g(0 * )= 0.
Простейшим градиентным методом является метод наискорейшего спуска, для которого формула (3.3) принимает вид:
0 ..+1 = 0.. - А. g...
Данный метод осуществляет поиск слишком медленно, поскольку градиент является локальной характеристикой и спуск, основанный на нем, ведет к частой смене направления (Банди, 1988). Более эффективны методы, использующие свойства квадратичных функций и корректирующие направление поиска по вторым производным минимизируемой функции. Базовым среди них является метод Ньютона с формулой поиска:
0 .+1 = 0. -[Н. I-1 g.,
где Н. — матрица Гессе (матрица вторых производных).
Комбинацию первых двух методов представляет собой метод Ньютона-Рафсона:
0.+1 = 0. -А. [Н. I-1 g.. (3.5)
Если начальное значение 0 0 задано достаточно близко к точке минимума 0 , то сходимость в данном методе будет быстрой, так как в локальной окрестности
154
этой точки квадратичная аппроксимация минимизируемой функции будет в общем случае хорошей.
Существуют специальные градиентные методы, относящиеся к классу нелинейных МНК (Демиденко, 1981), предназначенные для поиска с целевой функцией в виде суммы квадратов отклонений (СКО):
* (0)=! V,2 (в)=2[ - ^ (х; е)]2 = [ - г(е)Г [ - г(е)].
г г
В предположении, что величина ^ мала для всех t, матрица Гессе может быть аппроксимирована выражением:
Н5 - Н = 2РТР,
где Р = У е Г; Г = }={ (х,; б)}= Г(х; в); Уе — оператор градиента.
Градиент СКО также имеет специфический вид:
g, =У е 5 = -2РТ [у - р].
На основе указанной аппроксимации построены методы Ньютона-Гаусса и Левенберга-Марквардта (Демиденко, 1981). Первый метод производит поиск в соответствии с формулой:
Ö !+1 = 0 г - [HI-1 (g^) = 0 г + PP j"1 Pf |y - F ], (3.6)
где F = F(x; 0 J; p = Р(0г).
Главная проблема рассмотренного метода состоит в частом вырождении
T I T г 1
матрицы Рг. Рг., что не позволяет определить обратную к ней матрицу [р pj . Для устранения данной проблемы в указанную матрицу вводится поправочный член, что и составляет особенность метода Левенберга-Марквардта. В этом методе поиск ведется по формуле:
0 i+1 = 0i + |р + A.Mi j PT | - Fi ]. (3.7)
Поправка Левенберга является единичной матрицей M i = I, а поправка Мар-квардта — диагональной вида M i = diag(pT р). Масштабирующий множитель
pi > 0 обновляется специальным образом (Бокс, Дженкинс, 1974; Демиденко, 1981).
К сожалению, методы Ньютона-Гаусса и Левенберга-Марквардта пригодны только для стационарных гауссовских процессов. В противном случае аппроксимация матрицы Гессе (см. выше) становится очень неточной, и алгоритм не сходится. Таким образом, нелинейный МНК для целевой функции (3.1) может не подойти. Поэтому при решении задачи (3.2) следует также использовать и алгоритм для нелинейных функций общего вида.
Таким алгоритмом является метод Давидона-Флетчера-Пауэлла (ДФП), который также основан на квадратичной интерполяции, но не использует аппроксимацию матрицы Гессе. Более того, в данном методе обращение указанной матрицы вообще не требуется, так как матрица [H г j-1 рекуррентно аппроксимируется матрицей специального вида. При выполнении ряда не очень жестких ограничений в случае функции m переменных существует предел limGi = [Hm j-1 (Бан-
i^m
ди, 1988). Это первая из двух важных модификаций метода Ньютона. Вторая модификация связана с ведением поиска по сопряженным направлениям, которые являются наилучшими для поиска минимума квадратичных функций, и гаран-
155
тируют его нахождение для функции т переменных не более чем за т шагов (Банди, 1988). Основная формула поиска метода ДФП имеет вид:
0,.+1 = 0. -А.С.gг, (3.8)
где Gi — положительно определенная симметрическая матрица, обновляемая на каждом шаге итеративного цикла по формуле:
С !+1 = С. + А. + В
где матрицы А. и В. определены формулами:
А = у''в =-сс, и = g -g V =-1аР
А г Т , Т^ " г g г+1 g г Уг Л^г&г
V г и г и, Си г
В качестве начального значения G0 можно взять любую положительно определенную симметрическую матрицу, например единичную матрицу I. Хотя в целом метод ДФП осуществляет поиск медленнее, чем метод Ньютона-Рафсона, его сходимость к локальному минимуму имеет больше шансов.
4. Построение ОФК 1-го порядка для параметрической идентификации
Традиционный подход к определению градиента и матрицы Гессе для целевой функции — численный. Он основан на представлении первых и вторых производных в конечных разностях. Однако объем вычислений и время ожидания результата при таком прямом подходе непомерно большие, поскольку расчеты приходится делать по каждой компоненте вектора параметров 0, для каждого момента времени t и на каждом шаге итерационного цикла г. Поэтому целесообразно попытаться выразить искомые производные в аналитическом виде. Дифференцирование целевой функции (3.1) дает следующие формулы для градиента и матрицы Гессе:
Р1 = Vв1г = 2{[[ -V2.]еД,г + V^2г}/д2г;
г=о
н^ ^ т =(v е V е т ), ^ б е /,. =
= ,2 О, (б, О,,, + б, ^)- О,,, (б, О,,, + V, О,,,V„Р],)+ 2V, О,,, (^V,Р], + О,,V, у)")
2 п3
где I, =1 (Уо...„|0г); V . (0г) и А,,- = д (ег); бе = ^veТ) — матричный оператор вторых производных.
Далее, очевидно, необходимо найти аналитическое представление градиентов VвDti и Ve . Наиболее простым является представление искомых величин с помощью рекуррентных формул. Для их построения рекомендуется последовательно продифференцировать уравнения фильтра, выраженные в форме обновления (Справочник по прикладной статистике, 1990). Общая нелинейная модель
пространства состояний для одношагового предиктора Х+ в форме обновления легко выводится из исходной модели (1.8), (1.9) и уравнений фильтра, показанных на рис. 2:
х ++1 = f (Х ++ К (~ (; 0)+ и (. у ( = +; 0)+ ~ (
Отсюда непосредственно можно найти ковариацию обновляющего процесса:
Р = Н [^р-Х/+Q 4- ]н;г+R
где Н' = VХЬ, и Г = У.
Рассмотрим аналогичную модель для алгоритма (2.5)-(2.9)*:
) + — Т7+ -1-
пн-1 = 1 : иг; (4 1)
) + , У = пг + ^г
•г =
/ 'Л2
Р- + ^Р + (4.2)
1: 1 г-1
4-1' и р ' иш-
Для функции Риккера использовали значения (см. формулу (2.3)):
К =1 (( + кл; а ь)=[ ехР(а - Ьх, =я:+А>,; (43)
= [(;а,Ь)/дх,]== [( - Ьх,)ехр(а - Ьх,)]=(4.4)
Предварительно из формул (2.6), (2.7), (2.9) и (4.2) получили два полезных для дальнейших преобразований равенства:
р п 2 Кг = -р- и Кг = 1 - . (4.5)
г 2
•г
а . •
С помощью формул (4.1), (4.5) и правила дифференцирования сложной функции вывели требуемые рекуррентные соотношения для величин р {, • и их градиентов:
Р г :1 = Уг:1 - Пг+1 = Уг:1 - = Уг:1 -1« + Кг Рг;Ь) =
= Уг:1 -1 (У г - Р г + Кг Р г;Ь) = Уг :1 -1 (Уг -(1 - Кг )Р г;Ь) = (4.б) = Уг:1 - 1~г;
V у^ = [(1 - К, )Уу, -V, Уе К, ]; УеУ,2 = У^У,; (4.7)
Dt = а 2/(1 - *г); (4.8)
Уе Dt =[1/(1 - Кг )]Уе а2 +[а 2/(1 - К )2 Уе Кг; (4.9)
% = % (-(1 - К, )р,;аь )=[х, ехР(а - Ьх, )]х,=у,-(1-к, )р,; (410)
1~/=[ (х,;аь)дх,]х,=у,-(1-к,)р, =[(1 -Ьх,)ехР(а -Ьх,)]=у,-(1-к,)р,. (411)
Чтобы замкнуть полученные уравнения, необходимо получить рекуррентные соотношения для величины К и ее градиента по параметрам, а также задать все начальные значения и константы. Используя формулы (2.3), (2.6), (2.7) и (4.5), получили:
к = 12 Кг +' I ; ( )
Кг:1 = 2 ((2 г :1): 2 ; (4.12)
* Аргумент 0 для краткости опущен.
а 4 [2К V К 1+ а 2 V а 2 - а 2 V а 2
уу ^ аш I2 К г1 гу е1 уе К t Р^ °ту е ар а ру е ат
V е = [а 2 (2 К г +1)+ а 21 ; (4Л3)
V е р; = Ve К, -(1 - Кг) е V, ]+д Ve а + д е Ь; (4.14)
р'=[д2р(xt;а2^=л-(l-кt) = [(2 - ьх )ехР(а -ьх =л-(l-кtк; (4Л5)
д а 1 = 1; Э ЬР; =-& + =у-(-К, V . (4.1б)
Константные градиенты могут быть вычислены непосредственно:
V,, а = (1,0,0Д0)Т; V 0Ь = (0,1,0Д0)Т; V ао2р = 2^ (0,0,0,1,0)Т; V ^ = т (0,0,0,0,1)Т.
Остается задать начальные значения для величин р0 , Ve р0 , К0 и V е К0, и рекуррентный алгоритм (4.6)-(4.16) будет замкнут. В соответствии с начальным условием п0 = у0, предположили, что v0 = 0. Тогда Ve V 0 =(0,0,0,0,0 У. Из формулы (4.5) следует, что К0 = Р0/д2 и, следовательно,
Ve К 0 =
' 1 2Р ЛТ
0,0, ^г,- 2Р30,0
0
3
а т
Таким образом, построен алгоритм, который позволяет производить поиск минимума целевой функции (3.1) в удобной рекурсивной манере, вообще не прибегая к определению переменных состояния и дифференцированию.
5. Применение ОФК к прогнозированию подходов горбуши
Исходными данными для моделирования послужили нерестовые возвраты горбуши поколений четных лет в период 1970-2002 гг. на юго-восточное побережье о. Сахалин, включая зал. Терпения. Величина подхода в 2004 г. была определена приблизительно как сумма предварительных оценок официального вылова (7 млн экз.) и требуемого заполнения нерестилищ (7-11 млн экз.).
Прежде чем представить результаты параметрической настройки и применения ОФК в отношении указанных данных, следует сказать о внешнем сигнале ut в уравнении (2.1). Ввод переменной и обусловлен присутствием квазигармонического тренда в отклонениях фактических ежегодных возвратов горбуши от соответствующей кривой Риккера, в том числе для рассматриваемого района (Михеев, 2004а). Коэффициент корреляции между этим трендом и исходным рядом составил +0,5. Параметры тренда показаны в табл. 1. Внешний управляющий сигнал трактовали как результат влияния на системную динамику периодических факторов среды разного масштаба. Заметим, что в литературе приводятся убедительные примеры тесной связи между долгопериодными циклическими трендами в динамике обилия горбуши и рядом глобальных климатических и геофизических характеристик (Кляшторин, 2000).
В качестве алгоритма поиска оптимальных параметров для ОФК были последовательно использованы четыре метода из параграфа 3. Первоначально применили метод Ньютона-Гаусса (формула (3.6)), в соответствии с имеющимися в литературе рекомендациями (Справочник по прикладной статистике, 1990), однако вычисления часто либо не сходились к определенному результату, либо приводили к вырождению аппроксимирующей матрицы Гессе. В итоге при проведении расчетов этим методом потребовалось перебрать большое число пробных начальных значений. Учитывая, что дисперсия обновляющего процесса зависит от
времени (см. формулу (4.2)), предположили, что причина неустойчивости заключается в его нестационарности, поэтому попробовали алгоритм Ньютона-Рафсо-на (формула (3.5)). Однако этот алгоритм, хотя и сходился быстро в случае удачного пробного начального значения, регулярно приводил к вырождению матрицы Гессе. После того как установили, что обновляющий процесс стабилизируется, использовали метод Левенберга-Марквардта (формула (3.7)) в той модификации, как описано в работе Дж.Бокса и Г.Дженкинса (1974). Этот метод показал приемлемую сходимость и быстроту, однако в ряде случаев процесс вычислений прерывался, так как поправка в матрицу Гессе возрастала неограниченно. Поэтому, когда рассматриваемый метод не приводил к результату, поиск осуществляли с помощью метода ДФП (формула (3.8)). Хотя расчеты, производимые с использованием последнего метода, шли медленнее, они всегда завершались по критерию остановки.
Стартовые значения параметров модели (см. формулу (2.10)) задали следующим образом. Для параметров а0 и Ь0 пробные оценки получили с помощью
МНК для линеаризованной функции Риккера 1п(уг:1/уг )= а0 - Ь0уг. Величины ор и от перебирали последовательно с шагом 0,05, начиная от 0,1. Значение Р0 извлекали случайным способом с равной вероятностью из 10 отрезков, разбитых
в интервале от 0 до д^. Для начальных значений, приведенных в табл. 2, поиск
по методу ДФП завершился за 17 итераций. В табл. 2 показаны полученные оптимальные значения параметров.
Таблица 2
Оценки параметров для модели учета горбуши четных поколений на юго-востоке о. Сахалин (по результатам применения ОФК)
Table 2
Parameter estimations for account model of the even-year southeastern Sakhalin pink salmon generations (from the EFK application results)
„ Параметры фильтра
Показатель , r n
a b P0 a a
U p m
Начальные значения 0,897 0,057 0,038 0,200 0,200 Оценки ОФК 0,793 0,054 0,040 0,202 0,196
Как упоминалось выше, сходимость ОФК и выбор метода параметрической идентификации во многом зависят от стационарности обновляющего процесса. В связи с этим перед применением построенного алгоритма для прогнозирования целесообразно исследовать устойчивость дисперсии обновляющего процесса Dt (Справочник по прикладной статистике, 1990). Формула (4.2) показывает,
что поведение указанной величины определяется членом Ft'2 Pt. Из формул (2.6), (2.7) и (2.9) для ковариации ошибки вывели следующее рекуррентное соотношение:
2 т—т/2 т-) . 2 2
р = am Ft Pt + (51)
Pt+1 = n 2 , 2 . (51)
Ft Pt + + am
Пусть формулой (5.1) определена функция р+1 = f (р ). Нетрудно убедиться,
что данная функция обладает характеристиками: f (0)= а pö^/(^ 2 +
lim f (x)= а2 > f (0), f (x)> 0, f'(x)< 0, xe R0+ = [0,+~), т.е. является дифферен-
159
цируемои, монотонно возрастающей, всюду выпуклой и ограниченной на положительном луче вещественной прямоИ. Можно показать, что при любом фиксированном значении F' существует неподвижная точка отображения /: ^ , задаваемая формулоИ:
2 2
1 a + J1 а + ß; а = ^
2 . 2 % +
-г/2
; ß
-г/2
где Д*' — фиксированное значение . Как известно, у отображения F: ^ ^0+, определенного функциеИ Риккера F(х; а,Ь) = хехр(а -Ьх), для всех а е (0,2) существует неподвижная точка х = а/Ь с производноИ F»' = 1 — а (Скалецкая и др., 1979). Следовательно, для ОФК с функциеИ Риккера, у котороИ а < 2, дисперсия обновляющего процесса будет асимптотически стремиться к постоянному значению: ^ ^ ^ . В частном случае, когда а = 1, линия замещения функции Риккера проходит через ее максимум, и рассматриваемые величины принимают значения:
22
F = 0; P * =4^; D*
* p + ^m
22 ■^p +•
Применив построенныИ алгоритм ОФК к исходным данным, получили рекуррентную последовательность значениИ величины р. Эта последовательность показана на рис. 3. Отметим, что результаты расчетов полностью согласуются с результатами анализа и позволяют сделать вывод о стационарности обновляющего процесса. Кроме того, параметр а функции Риккера для рассматриваемых данных достаточно близок к единице (табл. 2).
II ПП'
и пгщ
«
ю
и U.04Ü
а П 042
о
к ■ I inn
и пп?
ц
а и 1124
и
р U.Olü
а
м U.012
о о.оио
6 S 10 12
Шаг итерации
Рис. 3. Поведение ковариации ошибки при применении ОФК к данным учета горбуши поколений четных лет на юго-востоке Сахалина
Fig. 3. Error covariation behaviour when applying EFK for account data of the even-year southeastern Sakhalin pink salmon generations
Построенная с помощью ОФК ненаблюдаемая динамика популяции для исходных данных и исходный ряд учтенной численности и последовательность од-ношаговых прогнозов с соответствующим доверительным интервалом показаны на рис. 4. Можно заметить, как по мере накопления данных происходит настройка фильтра, которая выражается во все лучшем "угадывании" тенденций в динамике подходов. Наряду с этим уменьшается, а затем и стабилизируется ко-вариация ошибки, характеристикой которой является доверительный интервал оценивания.
На рис. 5 показаны исходные и "отфильтрованные" данные и соответствующие кривые функции Риккера. Можно видеть почти полное совпадение этих
кривых в силу того, что применение ММП на основе ОФК мало повлияло на оценки параметров риккеровской функции, предварительно найденные с помощью МНК (табл. 2). Установленная близость оценок является следствием сходимости обновляющего процесса к стационарному состоянию и стабилизации фильтра (см. параграф 3). Тем не менее имеется заметное различие между точками на графике, отображающими исходные и отфильтрованные данные. Удаление шума из измерений привело к тому, что соответствующие точки легли ближе к теоретической кривой. Оставшиеся отклонения обусловлены вкладом известного внешнего сигнала.
Рис. 4. Результаты применения ОФК к данным учета горбуши поколений четных лет на юго-востоке Сахалина за период 1970-2004 гг.
Fig. 4. Results of the EFK application to account data of the even-year southeastern Sakhalin pink salmon generations in 1970-2004
д
о
г
в
а
д
о .з
х к
д э
о н
п л
ь
т
с
о 4-
н £
н
е
л
с
и
Ч
60-
40-"
20--
□ □
4»
□
— факт
— ФК
— МНК — ФК
О 20 40 60
Численность подхода в год t, млн экз.
Рис. 5. Кривые Риккера для горбуши поколений четных лет на юго-востоке Сахалина за период 1970-2002 гг.
Fig. 5. The Ricker's curves for the even-year southeastern Sakhalin pink salmon generations in 1970-2004
6. Обсуждение некоторых особенностей применения ФК
В ФК консолидированы многие ранее предложенные процедуры оценивания и реализован метод рекурсивного байесовского вывода для гауссовских рас-
161
пределений. Подход Калмана (1960) позволил построить алгоритм, обновляющий априорную информацию о системе по мере поступления данных и одновременно минимизирующий СКО. Рекурсивность стала важнейшим достоинством рассматриваемого фильтра и помогла решить ряд проблем, существующих в приложениях теории фильтрации до момента его создания. В частности, стало возможным отказаться от фильтров с конечной "памятью", что значительно повысило точность получаемых оценок (Kalman, 1960). Кроме того, в ФК проявились адаптивные свойства, выраженные в возможности априорной оценки точности получаемых результатов средствами самого алгоритма и его автонастройки в процессе обработки данных (Шильман, 1994; Eigbe et al., 1998). Наконец, благодаря своей рекурсивной структуре, ФК оказался незаменим для обработки данных в режиме "предиктор—корректор" в реальном времени. Еще одно новое качество, присущее ФК, — это представление уравнений фильтра в форме пространства состояний. Данное свойство расширило область приложения числовых фильтров со скалярных временных рядов на структурированные. В результате ФК стал пригодным для одновременной обработки всех доступных или представляющих интерес источников информации. Последнее, о чем следует упомянуть, — это способность ФК разделять источники шума и оценивать их характеристики — качество, весьма привлекательное при анализе промысловых данных.
В связи с применением ФК обычно указывается на ряд стандартных ограничений, а именно: динамика системы и измерений должна быть линейной, шум процесса и измерений является белым и гауссовским (Maybeck, 1979; Справочник по прикладной статистике, 1990). При соблюдении перечисленных выше условий оценка ФК является оптимальной (в смысле СКО) среди всех возможных оценок. Для обоснования применимости ФК в литературе (Maybeck, 1979) имеется некоторая аргументация в пользу того, что многие сложные системы, как природные, так и искусственного происхождения, являются почти линейными и содержат шумы, близкие по своим характеристикам к белому гауссовскому. Однако при ближнем рассмотрении требования, предъявляемые к ФК, либо не принципиальны, либо могут быть существенно снижены. Например, в случае негауссовско-го шума ФК имеет наименьшую ковариацию ошибки в классе линейных несмещенных фильтров (Maybeck, 1979). Заметим, что для вывода алгоритма на рис. 2 никаких предположений о форме плотности не требуется.
Общий подход к построению аппроксимирующих алгоритмов ФК для обработки негауссовских переменных и "цветного" шума изложен А.В.Балакришна-ном (1988). Уместен пример ОФК с коррелированным шумом, который в виде марковского процесса "случайного блуждания" вводился в параметр a функции Риккера при изучении влияния среды на репродуктивный потенциал популяций лососей (Peterman et al., 2003). В случае негауссовских плотностей лучшее качество в принципе обеспечивает ОФК, но его, как правило, гораздо труднее реализовать (Справочник по прикладной статистике, 1990). Для слабо нелинейных систем ОФК является асимптотически оптимальным. Вообще отказ от линейности системы не выводит ОФК из класса наилучших линейных оценок, но не гарантирует глобальной оптимальности в смысле достижения границы Крамера-Рао (Van der Merwe, 2004). Как было показано в параграфе 2, при настройке модели кардинальных методических отличий, связанных с нелинейностью, не возникает. Тем не менее в случае нелинейной модели всегда следует помнить о том, что алгоритм имеет дело с аппроксимацией. Поскольку при линеаризации уравнений фильтра отбрасываются все производные, начиная со второго порядка, применение ОФК 1-го порядка к сильно нелинейной модели часто ведет к смещению ковариации ошибки и ее быстрому росту. Тот же эффект при аппроксимации возникает в связи с неточностью в оценках параметров и начальных условиях. В свою очередь, смещение в оценке ковариационной матрицы может привести к расходимости фильтра. Ситуация обостряется, когда шум процесса является ма-
лым и в динамике превалирует нелинейность. Поэтому при работе с ОФК следует убедиться в стабилизации значений ковариации, как это было сделано в параграфе 5 (см. рис. 3). Если она отсутствует, то рекомендуется использовать фильтр с взвешенной памятью (Справочник по прикладной статистике, 1990). Предлагаются также адаптивные ФК, которые устойчивы к возмущениям начальных данных (Шильман, 1994). В последние годы появился целый класс точечных ФК (ТФК), которые справляются с проблемой расходимости в нелинейной модели с помощью специальных алгоритмов (Wan, Van der Merwe, 2000; Arulampalam et al., 2002; Van der Merwe, 2004). ТФК не зависят от формы плотностей и не нуждаются в разложении Тейлора. Достигается это с помощью стохастической аппроксимации и преобразований не самих переменных, а их плотностей (Julier, Uhlmann, 1997, 2004).
В зависимости от аппроксимации плотностей и способа вычисления приведенных выше интегралов (1.3) и (1.4) можно получить различные рекурсивные фильтры (Van der Merwe, 2004). Хотя рекурсивный байесовский вероятностный вывод прямо не использовался при построении оригинального ФК, для гауссовс-кой аппроксимации уравнения фильтра выводятся непосредственно из формулы Байеса (1.3) (Van der Merwe, 2004). Как известно, существуют различные по форме эквивалентные алгоритмы ФК. В связи с этим отметим, что автором получены два таких алгоритма ОФК для обновляющего процесса. Первый был представлен ранее (Михеев, 2003, 2004а), а второй — более простого вида — содержится в формулах (4.6)-(4.9). Эквивалентность указанных алгоритмов выражается в их взаимозаменяемости и сводимости. Кстати говоря, и замыкающие уравнения (4.10)-(4.16) у них одни и те же.
Полученный в данной работе ОФК не предназначен для параметрической идентификации в режиме on-line, поскольку параметры в нем рекурсивно не определяются. Так что в указанном смысле этот фильтр не является адаптивным (Шильман, 1994). Однако своих известных адаптивных свойств в отношении характеристик обновляющего процесса он, естественно, не утратил. Следует подчеркнуть, что при анализе промысловых данных второе качество, на наш взгляд, гораздо важнее, поскольку в большинстве модели рыболовства предполагают фиксированные параметры, по крайней мере на некоторый период времени. Поэтому и не ставилась цель построить алгоритм с пошаговой настройкой параметров, который был бы заметно сложнее в реализации. Если при поступлении наблюдений оценки параметров обновляются, они должны определяться, как и переменные состояния, с учетом только предыдущих значений и поступившего измерения. Решается такая задача методами совместного и сопряженного (дуального) оценивания. В первом случае исходная переменная расширяет число своих состояний за счет параметров, а уравнения модели остаются без изменения. Однако это приводит к тому, что даже линейная система становится нелинейной и требуется применять ОФК или ТФК (Van der Merwe, 2004). Во втором случае, как, например, у Петермана с соавторами (Peterman et al., 2003), для вектора параметров строится своя модель пространства состояний. В итоге образуются два ФК, которые работают последовательно: на вход первого подаются известные параметры и по ним оцениваются состояния; затем на вход второго подаются оцененные состояния, и по ним оцениваются параметры; цикл повторяется.
Какой бы ни была структура ФК, реализация его алгоритма требует решения соответствующей оптимизационной задачи. И не важно, в каком контексте, байесовском или частотном, она сформулирована. Традиционные градиентные методы, использованные в настоящей работе, были подробно рассмотрены в параграфе 2. В качестве примера их приложения к оцениванию параметров ФК можно назвать применение метода квази-Ньютона ДФП в исследовании стохастических трендов для различных моделей ресурсопользования (Perrings, Stern, 2000). Добавим, что, по мнению некоторых авторов (Meyer, Millar, 1999; Millar, Meyer,
163
2000), необходимость в градиентных методах и ММП является одним из главных недостатков ФК. Частично с этим можно согласиться, поскольку проблема действительно есть. Во многом она связана с формой функции логарифмического правдоподобия, по поводу которой существуют противоположные мнения. Например, что она "плоская" (Schnute, 1994) или что она, как правило, многомодальная, может иметь разрывы на границе и "седловые" точки (Sekhon, Mebane, 1998). В последнем случае при применении градиентных методов результаты будут сильно зависеть от выбора начальных значений для искомых параметров. Как обсуждалось выше, даже в случае такой простой модели, как риккеровская, это действительно случается. Тем не менее приведенное выше утверждение Майера и Миллара (Meyer, Millar, 1999) о непригодности ФК к ряду рыбопромысловых моделей в связи с проблемами оптимизации параметров является, на наш взгляд, чрезмерно сильным. Теория фильтрации постоянно развивается, и, как неоднократно упоминалось, многие из стоявших ранее проблем утратили свою актуальность.
Получили развитие и методы оптимизационного поиска. Далее вкратце обозначим другие появившиеся подходы к решению данного вопроса. К настоящему времени сформировалось еще два принципиально разных направления: методы, базирующиеся на автоматическом дифференцировании (Griewank, Corliss, 1991), и методы вероятностного эволюционного поиска (генетические алгоритмы) (Воро-новский и др., 1997). Сегодня названные методы оцениваются, в целом, как более эффективные, чем градиентные (Schnute et al., 1998; Sekhon, Mebane, 1998). Однако сходимость многих генетических алгоритмов строго не доказана, и нахождение глобального экстремума с их помощью не гарантируется, хотя шансы значительно выше. С другой стороны, градиентные методы работают значительно быстрее, и в случае малого числа параметров их применение часто является достаточным. Отмеченное обстоятельство стало поводом для создания методов гибридного поиска, сочетающих в себе локальную быстроту градиентных и глобальную эффективность эволюционных алгоритмов (Вороновский и др., 1997). Методы поиска, основанные на автоматическом дифференцировании и генетических алгоритмах, внедряются и в рыбопромысловые модели (Fournier, 1994; Saila, 1996; Chen et al., 2000, 2002).
В отношении предсказательной способности ОФК на основе предложенной модели учета отметим следующее. Во-первых, несмотря на значительную неопределенность в предварительной оценке фактического подхода в 2004 г., тенденция была угадана точно, как, впрочем, и для десяти предыдущих лет (см. рис.4). Да и возможная погрешность колеблется в допустимых пределах: от +3 до +30 % при оценке ОФК в 18,5 млн экз. Во-вторых, прогноз для генераций нечетных лет этого же стада на 2003 г. имел ошибку всего в +2,9 % (Михеев, 2003, 2004а).
В настоящей работе представлен новый алгоритм оптимизации параметров для ОФК 1-го порядка, базирующийся на аналитическом вычислении градиентов логарифмического правдоподобия. Проведенные исследования позволили обосновать сходимость указанного фильтра для модели учета одновозрастной популяции, описываемой функцией Риккера. Полученные результаты анализа были подтверждены при моделировании динамики нерестовых возвратов горбуши юго-восточного Сахалина за период 1970-2002 гг. С помощью построенного ОФК реконструировали динамику подходов названного промыслового стада и дали ее прогноз на 2004 г.
Согласованность результатов анализа и моделирования в отношении сходимости представленного ОФК позволяет предложить его в качестве одного из методов прогнозирования возвратов как горбуши, так и других промысловых видов. Другая перспектива для этого метода открывается в направлении исследований влияния внешней среды на динамику численности тихоокеанских лососей в долгопериодном и краткосрочном аспектах. Следует еще раз сказать, что боль-
шинство ограничений ОФК, относящихся к типу шума и линейности системы, сегодня уже сняты или существенно ослаблены. Большинство проблем, встречающихся в практике оптимизации оценок ОФК, к настоящему времени преодолены, оставшиеся потенциально решаемы. С другой стороны, рассматриваемый фильтр имеет ряд принципиальных позитивных свойств, среди которых разделение и оценка характеристик источников шума, оптимальная фильтрация полезного сигнала из данных наблюдений, адаптация к поступающей информации. Они делают данный метод эффективным средством обработки некачественной промысловой статистики, что нельзя не учитывать при определении текущих и перспективных величин запаса. Прикладная значимость рассмотренного метода прогнозирования состоит также и в том, что он позволяет постепенно оценить ошибки учета, что особенно актуально при масштабном браконьерстве.
Литература
Балакришнан A.B. Теория фильтрации Калмана. — М.: Мир, 1988. — 168 с.
Банди Б. Методы оптимизации. Вводный курс. — М.: Радио и связь, 1988. — 128 с.
Бокс Дж., Дженкинс Г. Анализ временных рядов. Прогноз и управление. — М.: Мир, 1974. — Вып. 1. — 406 с.
Брандт 3. Анализ данных. Статистические и вычислительные методы для научных работников и инженеров. — М.: Мир; ООО "Издательство АСТ", 2003. — 686 с.
Вероятность и математическая статистика: Энциклопедия. — М.: Бол. Рос. Энциклопедия, 1999. — 910 с.
Волков И.К., Зуев С.М., Цветкова Г.М. Случайные процессы. — М.: Изд-во МГТУ им. Н.Э.Баумана, 2003. — 448 с.
Вороновский Г.К., Махотило К.В., Петрашев С.Н., Сергеев С.Н. Генетические алгоритмы, искусственные нейронные сети и проблемы виртуальной реальности. — Харьков: Основа, 1997. — 112 с.
Демиденко Е.З. Линейная и нелинейная регрессия. — М.: Финансы и статистика, 1981. — 302 с.
Кляшторин Л.Б. Тихоокеанские лососи: климат и динамика запасов // Рыб. хоз-во. — 2000. — № 4. — С. 32-34.
Колмогоров А.Н. Интерполирование и экстраполирование стационарных случайных последовательностей // Изв. АН СССР. Сер. Мат. — 1941. — Т. 5.
Корельский В.Ф., Макоедов А.Н. Регулирование конкурентных отношений при промысле водных биологических ресурсов // Вопр. рыб-ва. — 2004. — Т. 5, № 3. — С. 385-393.
Крамер Г. Математические методы статистики. — М.: Мир, 1975. — 648 с.
Математический энциклопедический словарь / Гл. ред. Ю.В.Прохоров. — М.: Сов. энциклопедия, 1988. — 847 с.
Михеев A.A. Метод среднесрочного прогнозирования нерестовых подходов горбуши в заливе Анива и на восточном побережье о. Сахалин: Отчет по гранту Администрации Сахалинской обл., проект № 2-2.1-01. — Южно-Сахалинск: СахНИРО, 2003. — 61 с.
Михеев A.A. Моделирование стохастических процессов в эксплуатируемых популяциях рыб и беспозвоночных (на примере горбуши Oncorhynchus gorbuscha и синего краба Paralithodes platypus восточного шельфа Сахалина): Дис. ... канд. биол. наук. — Южно-Сахалинск: СахНИРО, 2004а. — 192 с.
Михеев A.A. Применение фильтра Калмана к модели Риккера "запас-пополнение" // Тез. докл. семинара "Математическое моделирование и информационные технологии в исследованиях биоресурсов Мирового океана". — Владивосток: ТИНРО-центр, 2004б. — С. 29-32.
Скалецкая Е.И., Фрисман Е.Я., Шапиро A.n. Дискретные модели динамики численности популяций и оптимизации промысла. — М.: Наука, 1979. — 166 с.
Справочник по прикладной статистике: В 2 т. — М.: Финансы и статистика,
1990. — Т. 2. — 526 с.
Шильман С.В. Адаптивные фильтры Калмана // ДАН СССР. — 1994. — Т. 338, № 6. — С. 742-744.
Arulampalam M.S., Maskell S., Gordon N., Clapp T. A tutorial on particle filters for on-line nonlinear // Non-Gaussian Bayesian tracking. IEEE Transactions on Signal Processing. — 2002. — Vol. 50, № 2. — P. 174-188.
Chen D.G., Hargreaves N.B., Ware D.M., Liu Y. A fuzzy logic model with genetic algorithm for analyzing fish stock-recruitment relationships // Can. J. Fish. Aquat. Sci. — 2000. — Vol. 57. — P. 1878-1887.
Chen D.G., Irvine J.R., Cass A.J. Incorporating Allee effects in fish stock-recruitment models and applications for determining reference points // Can. J. Fish. Aquat. Sci. — 2002. — Vol. 59. — P. 242-249.
Chen Y., Wilson C. A simulation study to evaluate impacts of uncertainty on the assessment of American lobster fishery in the Gulf of Maine // Can. J. Fish. Aquat. Sci. — 2002. — Vol. 59. — P. 1394-1403.
Eigbe U., Beck M.B., Wheater H.S., Hirano F. Kalman filtering in groundwater flow modeling: problems and prospects // Stochastic Hydrology and Hydraulics. — 1998. — № 12. — P. 15-32.
Fournier D.A. An introduction to AD Model Builder, for use in nonlinear modelling and statistics. — Canada, Nanaimo: Otter Research Ltd., 1994. — 51 p.
Francis R.I.C.C., Shotton R. "Risk" in fisheries management: a review // Can. J. Fish. Aquat. Sci. — 1997. — Vol. 54. — P. 1699-1715.
Griewank A., Corliss C.F. (Ed.) Automatic differentiation of algorithms: theory, implementation, and application. — Philadelphia, Pa.: Society for Industrial and Applied Mathematics, 1991. — 353 p.
Gudmundsson G. Time series analysis of catch-at-age observations // Appl. Stat. — 1994. — Vol. 43, № 1. — P. 117-126.
Gudmundsson G. Joint time series analysis of catch-at-age and CPUE data // Fish. Stock Assess. Models.: Alaska Sea Grant College Program Report № AK-SG-98-01. — Fairbanks: Univ. of Alaska, 1998. — P. 199-218.
Gudmundsson G. Time series analysis of abundance indices of young fish. ICES // J. Mar. Sci. — 2004. — Vol. 61. — P. 176-183.
Hilborn R., Walters C.J. Quantitative fisheries stock assessment: choice, dynamics and uncertainty. — N.Y.: Chapman and Hall, 1992. — 570 p.
Julier S.J., Uhlmann J.K. A new extension of the Kalman filter to nonlinear systems // Proc. SPIE. — Orlando, FL, USA: Int. Soc. Opt. Eng., 1997. — Vol. 3068. — P. 182-193.
Julier S.J., Uhlmann J.K. Unscented filtering and nonlinear estimation // Proc. of the IEEE. — 2004. — Vol. 92, № 3. — P. 401-422.
Kalman R.E. A new approach to linear filtering and prediction problems // J. Basic Eng. — 1960. — Vol. 82. — P. 34-45.
Kalman R.E., Bucy R.S. New results in linear filtering and prediction theory // J. Basic Eng. — 1961. — Vol. 83. — P. 95-108.
Kimura D.K., Balsiger J.W., Ito D.H. Kalman filtering the delay-difference equation: practical approaches and simulations // Fish. Bull. — 1996. — Vol. 94. — P. 678-691.
Maybeck P.S. Stochastic models, estimation, and control. — N.Y.: Academic Press, 1979. — Vol. 1.
Meyer R., Millar R.B. BUGS in Bayesian stock assessments // Can. J. Fish. Aquat. Sci. — 1999. — Vol. 56. — P. 1078-1087.
Millar R.B., Meyer R. Bayesian state-space modeling of age-structured data: fitting a model is just the beginning // Can. J. Fish. Aquat. Sci. — 2000. — Vol. 57. — P. 43-50.
Pella J.J. Utility of structural time series models and the Kalman filter for predicting consequences of fishery actions // Proc. Intern. Sympos. on Management strategies for exploited fish populations: Alaska Sea Grant College Program Report № AK-SG-93-02. — Fairbanks, Alaska: Univ. of Alaska, 1993. — P. 571-594.
Perrings C.A., Stern D.I. Modelling loss of resilience in agroecosystems: to rangelands degradation in Botswana // Environmental and Resource Economics. — 2000. — Vol. 16. — P. 185-210.
Peterman R.M., Pyper B.J., Grout J.A. Comparison of parameter estimation methods for detecting climate-induced changes in productivity of Pacific salmon (Oncorhynchus spp.) // Can. J. Fish. Aquat. Sci. — 2000. — Vol. 57. — P. 181-191.
Peterman R.M., Pyper B.J., MacGregor B.W. Use of the Kalman filter to reconstruct historical trends in productivity of Bristol Bay sockeye salmon (Oncorhynchus ner-ka) // Can. J. Fish. Aquat. Sci. — 2003. — Vol. 60. — P. 809-824.
166
Reed W.J., Simons C.M. Analyzing catch-effort data by means of the Kalman filter // Can. J. Fish. Aquat. Sci. — 1996. — Vol. 53. — P. 2157-2166.
Saila S.B. Guide to some computerized artificial intelligence methods // Computers in fisheries research / Ed. by B.A.Megrey and E.Moksness. — L.: Chapman and Hall, 1996. — P. 8-40.
Schnute J.T. A general framework for developing sequential fisheries models // Can. J. Fish. Aquat. Sci. — 1994. — Vol. 51. — P. 1676-1688.
Schnute J.T. The importance of noise in fish population models // Fish. Res. —
1991. — Vol. 11. — P. 197-123.
Schnute J.T., Kronlund A.R. Estimating salmon stock-recruitment relationships from catch and escapement data // Can. J. Fish. Aquat. Sci. — 2002. — Vol. 59. — P. 433-449.
Schnute J.T., Richards L.J. Use and abuse of fishery models // Can. J. Fish. Aquat. Sci. — 2001. — Vol. 58. — P. 1-8.
Schnute J.T., Richards L.J., Olsen N. Statistics, software, and fish stock assessment // Fish. Stock Assess. Models.: Alaska Sea Grant College Program Report № AK-SG-98-01. — Fairbanks: Univ. of Alaska, 1998. — P. 171-184.
Sekhon J.S., Mebane W.R. Genetic optimization using derivatives // Political Analysis. — 1998. — Vol. 7. — P. 187-210.
Sorenson H.W. Least-square estimation: from Gauss to Kalman // IEEE Spectrum. — 1970. — Vol. 7. — P. 63-68.
Sullivan P.J. A Kalman filter approach to catch-at-length analysis // Biometrics. —
1992. — Vol. 48. — P. 237-257.
Van der Merwe R. Sigma-point Kalman filters for probabilistic inference in dynamic state-space models. — D.Ph.Th., OGI School of Science & Engineering, Oregon Health & Science University, 2004. — 378 p.
Wan E.A., Van der Merwe R. The unscented Kalman filter for nonlinear estimation // Proc. IEEE Sympos. on Adaptive systems for signal processing communications and control (AS-SPCC). — Lake Louise, Alberta, Canada, 2000. — P. 153-158.
Welch G., Bishop G. An introduction to the Kalman filter. — ACM Inc., 2001. —
47 p.
Wiener N. The extrapolation, interpolation and smoothing stationary time series. — N.Y.: Wiley, 1949.
Поступила в редакцию 15.11.05 г.
