Верификация процесса конвективной диффузии на основе анализа многомерных временных рядов © М.Г. Матвеев © Е.А. Сирота © М.Е. Семенов © А.В. Копытин Воронежский государственный университет, Воронеж, Россия mgmatveev@yandex.ru atoris@list.ru Аннотация. Предложен метод верификации процесса конвективной диффузии на основе анализа многомерных временных рядов. Он основан на сравнительном анализе конечно-разностных представлений исследуемой модели и авторегрессионного описания временных рядов наблюдений. Метод включает получение МНК-оценок параметров многомерной авторегрессии и построение вариантов систем алгебраических уравнений, связывающих оценки авторегрессии и параметры соответствующих дифференциальных уравнений. Система алгебраических уравнений, которой удовлетворяют полученные оценки, определяет структуру модели и соответствующие значения параметров дифференциального уравнения. Приведен численный пример идентификации процесса изменения температуры атмосферного воздуха. Ключевые слова: уравнения в частных производных; структурная идентификация; параметрическая идентификация; многомерная авторегрессия; статистические гипотезы. Verification of the Convective Diffusion Process Based on the Analysis of Multidimensional Time Series © M.G. Matveev © E.A. Sirota © M.E. Semenov © A.V. Kopytin Voronezh State University, Voronezh, Russia mgmatveev@yandex.ru atoris@list.ru Abstract. A method for verifying the process of convective diffusion based on the analysis of multidimensional time series is proposed. The proposed method is based on a comparative analysis of finite- difference representations of the model under study and an autoregressive description of time series of observations. The method includes the LSM (Least square method) estimates of the parameters of multidimensional autoregression and the construction of versions of systems of algebraic equations connecting the estimates of autoregression and the parameters of the corresponding differential equations. The system of algebraic equations satisfied by the obtained estimates determines the structure of the model and the corresponding values of the parameters of differential equation. A numerical example of identifying the process of changing the temperature of atmospheric air is given. Keywords: partial differential equations, structural identification, parametric identification, multidimensional autoregression, statistical hypothesis. механизмов процессов могут служить методы 1 Введение анализа временных рядов, образующихся при Математические модели физических процессов наблюдении за параметрами функционирования этих основаны на априорном предположении о механизме процессов [2, 3]. Предлагаемый метод основан на функционирования этих процессов [1]. Однако такие сравнительном анализе конечно-разностных предположения не всегда в достаточной степени представлений исследуемой модели процесса и обоснованы. Ошибки спецификации модели авторегрессионнного описания временных рядов возникают как на структурном, так и на наблюдений. параметрическом уровнях. Эффективным Будем рассматривать широкий класс инструментом проверки предполагаемых пространственно-распределенных динамических систем, для которых характерны диффузионные Труды XIX Международной конференции процессы, процессы адвекции или их сочетание. «Аналитика и управление данными в областях с Соответствующее дифференциальное уравнение в интенсивным использованием данных» (DAMDID/ RCDL’2017), Москва, Россия, 10–13 октября 2017 года 354 частных производных с начальными и граничными 𝑖 = 2, … 𝑛 − 2; 𝑡 = 0,1, …, условиями имеет следующий общий вид: 0 𝑦 0 = (𝑦𝑖−1 ; 𝑦𝑖0 ; 𝑦𝑖+1 0 ), 𝜕𝑦 𝜕𝑦 𝜕2 𝑦 𝜕𝑡 + 𝜗 𝜕𝑙 = 𝐷 𝜕𝑙2 , (1) 0 𝑦𝑖−1 = (𝑦𝑖−1 1 ; 𝑦𝑖−1 𝑘 ; … ; 𝑦𝑖−1 ), 𝑦(0, 𝑙) = 𝜑(𝑙), 𝑦(𝑡, 𝑙min ) = 𝑓1 (𝑡), 𝑦(𝑡, 𝑙max ) = 𝑓1 (𝑡), 0 𝑦𝑖+1 = (𝑦𝑖+1 1 ; 𝑦𝑖+1 𝑘 ; … ; 𝑦𝑖+1 ), где 𝜗  скорость адвекции, 𝐷  коэффициент где 𝛼𝑖 – коэффициенты, вид которых определяется диффузии, 𝑙  пространственная координата. вариантом структуры процессов. Источником информации о поведении системы Начальные и граничные условия задаются как являются данные натурных измерений переменной результаты натурных измерений, поэтому 𝑦𝑖𝑡 с погрешностью 𝜀𝑖𝑡 в виде «белого шума» – 𝑥𝑖𝑡 = 0 𝑦 0 = 𝑥 0 = (𝑥𝑖−1 ; 𝑥𝑖0 ; 𝑥𝑖+1 0 ), 𝑦𝑖𝑡 + 𝜀𝑖𝑡 в последовательные моменты времени 𝑡 = 0 1 𝑘 0,1, …⁡в узлах одномерной пространственной 𝑦𝑖−1 = 𝑥𝑖−1 = (𝑥𝑖−1 ; 𝑥𝑖−1 ; … ; 𝑥𝑖−1 ), регулярной сетки⁡𝑖 = 0,1, … , 𝑛, т. е. многомерный 0 1 𝑘 𝑦𝑖+1 = 𝑥𝑖+1 = (𝑥𝑖+1 ; 𝑥𝑖+1 ; … ; 𝑥𝑖+1 ). временной ряд. Рассмотрение одномерной сетки ничем не ограничивает дальнейшие исследования, Заметим, что в соответствии со свойством зато позволяет избежать громоздких построений, консервативности [4] сумма коэффициентов правой характерных для плоских и объемных пространств. части выражений (2)(4) равна единице. Пусть в рассматриваемой системе могут протекать Соответственно ∑3𝑗=1 𝛼𝑗 = 1. процессы диффузии и адвекции. Утверждать, что поведение системы определяется одним из 3 Построение модели многомерного указанных процессов или действуют оба процесса временного ряда и ее сравнение с одновременно, нет достаточных оснований. Также разностной схемой неизвестны параметры (1), которые рассматриваются как константы, полученные в результате усреднения Моделирование многомерного временного ряда по времени. будем осуществлять в классе линейных Задача заключается в верификации процессов стохастических моделей авторегрессии [3]. конвективной диффузии на основе анализа Допустим, что в каждом узле регулярной сетки многомерных временных рядов и разработке протекает марковский процесс без последействия, и алгоритмов структурной и параметрической временные ряды в смежных узлах имеют высокие идентификации механистической модели с значения коэффициентов линейной корреляции, что постоянными коэффициентами по наблюдаемым обуславливает рассмотрение многомерных рядов. Такие допущения позволяют специфицировать значениям 𝑥𝑖𝑡 . стохастическую модель в i -м узле в виде 2 Разностные схемы для вариантов 𝑥𝑖𝑡+1 уравнений механистической модели 𝑥̃𝑖𝑡+1 = 𝑀 ( 𝑡 , 𝑥𝑖𝑡 , 𝑥𝑖+1 𝑡 𝑡 ) = 𝛼1 𝑥𝑖−1 + 𝛼2 𝑥𝑖𝑡 + 𝛼3 𝑥𝑖+1 𝑡 ; 𝑥𝑖−1 Для решения задачи составим явные трехточечные разностные схемы для уравнений 𝑖 = 2, … 𝑛 − 2; ⁡⁡⁡⁡⁡⁡𝑡 = 0,1, …, (6) каждого из вариантов структуры процессов: где⁡⁡𝛼1 , 𝛼2 , 𝛼3 – оценки параметров авторегресcии. диффузия и адвекция 𝑦𝑖𝑡+1 −𝑦𝑖𝑡 𝑦𝑡 −𝑦𝑡 𝑦𝑡 −2𝑦𝑡 +𝑦𝑖−1 𝑡 Вычисления по выражениям (5) и (6) существенно ∆𝑡 + 𝜗 𝑖+12∆𝑙 𝑖−1 = 𝐷 𝑖+1 ∆𝑙2𝑖 , (2) различаются, несмотря на их кажущееся сходство. Выражение (5) представляет собой рекуррентную 𝑦𝑖𝑡+1 = (𝑏1 + 𝑏2 )𝑦𝑖−1 𝑡 + (1 − 2𝑏2)𝑦𝑖𝑡 + (𝑏2 − 𝑏1 )𝑦𝑖+1 𝑡 ; формулу для вычисления модельных значений 𝜗∆𝑡 𝐷∆𝑡 𝑏1 = ;⁡⁡⁡⁡⁡𝑏2 = 2 ; переменной у в трех попарно смежных узлах сетки, 2∆𝑙 ∆𝑙 с условной устойчивостью, следовательно, и диффузия в неподвижной среде сходимостью к решению соответствующего 𝑦𝑖𝑡+1 − 𝑦𝑖𝑡 𝑡 𝑦𝑖+1 − 2𝑦𝑖𝑡 ⁡+ 𝑦𝑖−1 𝑡 =𝐷 ; дифференциального уравнения. Выражение (6) ∆𝑡 ∆𝑙2 позволяет вычислять модельные значения случайной переменной х в тех же узлах, но автономно для 𝑦𝑖𝑡+1 = 𝑏2 𝑦𝑖−1 𝑡 + (1 − 2𝑏2)𝑦𝑖𝑡 + 𝑏2 𝑦𝑖+1 𝑡 ; (3) каждого узла в каждый момент времени, как это адвекция определяется формулой (6). Различие в 𝑦𝑖𝑡+1 −𝑦𝑖𝑡 𝑦𝑡 −𝑦𝑡 вычислительных алгоритмах (5) и (6) затрудняет ∆𝑡 + 𝜗 𝑖+12∆𝑙 𝑖−1 = 0, (4) сравнение результатов вычислений. В частности, для 𝑦𝑖𝑡+1 = 𝑏1𝑦𝑖−1 𝑡 + 𝑦𝑖𝑡 − 𝑏1 𝑦𝑖+1 𝑡 . целей проводимого исследования важно показать, Представим разностные схемы (2)(4) в виде что математическое ожидание оценок а в обобщенной рекуррентной зависимости для авторегрессии (6) равно параметрам  в разностной произвольного узла i с начальными и граничными схеме (5). условиями: Найдем математические ожидания левой и правой части авторегрессионной зависимости, учитывая, что 𝑦𝑖𝑡+1 = 𝛼1 𝑦𝑖−1 𝑡 + 𝛼2 𝑦𝑖𝑡 + 𝛼3 𝑦𝑖+1 𝑡 (5) 355 математическое ожидание условного температурой атмосферного воздуха над заданным математического ожидания случайной величины участком земной поверхности. есть математическое ожидание этой случайной Для экспериментальной апробации величины: использовались статистические данные реанализа 𝑀(𝑥𝑖𝑡+1 ) = 𝑀(𝑎1 )𝑀(𝑥𝑖−1 𝑡 ) + 𝑀(𝑎2 )𝑀(𝑥𝑖𝑡 ) + параметров атмосферы за 2012 год [10]. Эти данные 𝑡 +⁡𝑀(𝑎3 )𝑀(𝑥𝑖+1 ). представляют собой ежедневные значения температуры в узлах плоской регулярной сетки с Полученное выражение можно переписать с учетом шагом 2,5°. Рассматривались результаты граничных условий и очевидного условия наблюдений температуры при геопотенциале 300 𝑀(𝑥𝑖𝑡+1 ) = 𝑀(𝑦𝑖𝑡+1 ) + 𝑀(𝜀𝑖𝑡+1 ) = 𝑦𝑖𝑡+1 ГПа в узле сетки с координатами 35° северной широты; 7,5° восточной долготы. По указанным в виде данным была построена одномерная модель 𝑦𝑖𝑡+1 = 𝑀(𝑎1 )𝑦𝑖−1 𝑡 + 𝑀(𝑎2 )𝑦𝑖𝑡 + 𝑀(𝑎3 )𝑦𝑖+1 𝑡 , множественной авторегрессии 𝑡+1 что доказывает равенство математических ожиданий 𝑥̃7,5 = 1,03𝑥5𝑡 − 1,14𝑥7,5 𝑡 𝑡 + 1,11𝑥10 𝑡 𝑥10 , (10) оценок параметров авторегрессии и параметров  где 𝑡 = 0, 1, …. Заметим, что сумма коэффициентов в разностной схемы (5). Полученный результат выражении (10) равна единице, что можно позволяет приравнять оценки параметров рассматривать как признак того, что временной ряд авторегрессии и параметров  и из полученных адекватно отражает динамику температуры, уравнений найти параметры 𝑏1 и 𝑏2. описанную разностной схемой (5). Для оценки Трем⁡ вариантом⁡ сочетания⁡ процессов⁡ будут⁡ процессов, оказывающих доминирующее влияние на соответствовать⁡ три⁡ варианта⁡ систем⁡ динамику температуры, проведем статистический алгебраических⁡ уравнений⁡ относительно⁡ анализ систем (7)–(9). параметров⁡𝑏1 и⁡𝑏2 разностных⁡схем⁡(2)–(4): Стандартные ошибки оценок равны 𝜎𝑎1 = для диффузии и адвекции 0,14;⁡⁡𝜎𝑎2 = 0,25; 𝜎𝑎3 = 0,13. Подставим 𝑏1 + 𝑏2 = 𝑎1 , полученные оценки 𝑎 в систему (7), что {1 − 2𝑏2 = 𝑎2 , (7) соответствует принятию гипотезы о совместном 𝑏2 − 𝑏1 = 𝑎3 ; влиянии на температурную динамику процессов для диффузии диффузии и адвекции. Допустим, что ошибки оценок 𝑏2 = 𝑎1 , 𝑎 равномерно распределены между коэффициентами {1 − 2𝑏2 = 𝑎2 , (8) 𝑏. Тогда последние два уравнения системы (7) 𝑏2 = 𝑎3 ; позволяют вычислить оценки 𝑏2 = 1,07 со для адвекции стандартной ошибкой 𝜎𝑏2 = 0,07 и⁡𝑏1 = −0,04 со 𝑏1 = 𝑎1 , стандартной ошибкой 𝜎𝑏1 = 0,125. Последнее {1 = 𝑎2 , (9) равенство означает, что 𝑏1 можно считать равным −𝑏1 = 𝑎3 . нулю. Поскольку коэффициент 𝑏1 отвечает за В полученных уравнениях в правой части всюду адвекцию (см. раздел 2), нулевая гипотеза о стоят оценки, т. е. значения случайных величин. совместном влиянии процессов диффузии и Следовательно, решение каждого из уравнений адвекции отклоняется. следует проводить как проверку соответствующей Подставим полученные оценки в систему (8), что статистической гипотезы в предположении соответствует принятию нулевой гипотезы о влиянии нормального распределения оценок. Вариант диффузии. В этом случае нулевая гипотеза может системы уравнений, согласующийся с быть принята, если подтвердятся гипотезы о соответствующей гипотезой, определит адекватный равенстве 𝑎1 и 𝑎2 и тождественности второго процесс. Решение с проверкой статистической равенства системы (8). Рассмотрим гипотезу о гипотезы показано в примере, приведенном в равенстве нулю разности 𝑎1 и 𝑎2 . Эта разность следующем разделе. составляет 0,08 при стандартной ошибке 𝑎1 и 𝜎𝑎2−𝑎1 = 0.13 + 0.14 = 0.27, что дает основание 4 Пример реализации предложенной принять гипотезу о равенстве и улучшить оценку методики путем вычисления среднего значения 𝑏2 = (1.03 + Диффузия и адвекция являются характерными 1.11)/2 = 1.07. Для проверки гипотезы о процессами, определяющими изменение состояния тождественности второго равенства системы (8) земной атмосферы. В частности, именно эти подставим в него полученные оценки 𝑎2 и 𝑏2 и процессы определяют динамику температурных получим –1,14=–1,14, что при стандартной ошибке полей на заданных изобарических поверхностях [8, 𝜎𝑎2 = 0.25 позволяет принять гипотезу о 9]. Актуальной задачей метеорологии является тождественности, следовательно, и нулевую определение структуры и параметров атмосферных гипотезу о влиянии на динамику температуры процессов. Покажем, что решение этой задачи процесса диффузии. возможно, по крайней мере, для среднегодовых Подстановка оценок 𝑎 в систему (9) для процесса значений, по результатам регулярных наблюдений за адвекции позволяет сразу отклонить гипотезу о 356 влиянии на динамику температуры процесса приближения производных в условиях помех [13]. адвекции из-за наличия противоречивых знаков Различные статистические методы были оценок. разработаны для оценки параметров моделей, Окончательный вывод: на динамику температуры описываемых обыкновенными дифференциальными атмосферы в рассматриваемой точке в среднем в уравнениями. Так, например, в [14–16] предложены течение года доминирующее влияние оказывали иерархические байесовские подходы к этой процессы атмосферной диффузии. проблеме. Эти методы требуют неоднократного численного решения соответствующих Чтобы убедиться в адекватности результатов обыкновенных дифференциальных уравнений, что в структурной и параметрической идентификации, ряде случаев требует значительного временного было найдено численное решение уравнения ресурса. Для оценки постоянных параметров модели диффузии по неявной разностной схеме, поскольку в в [17] предложен обобщенный сглаживающий данном случае условие устойчивости явной схемы не подход, основанный на идеях метода максимального выполняется. Уровни наблюдаемого временного правдоподобия. В работе [18] рассматривается ряда температур, результаты моделирования с каскадный метод оценки изменяющихся по времени использованием авторегрессии (10) и решение параметров. Этот метод оценивает параметры путем уравнения диффузии с использованием неявной оптимизации определенного критерия, однако разностной схемы с найденным значением параметра достижение глобального минимума проблематично с 𝑏2 = 2.02 представлены на рис. 1. вычислительной точки зрения. Другим подходом к оцениванию параметров обыкновенных дифференциальных уравнений (см. [19]) является двухэтапный метод, в котором на первом этапе, с использованием сглаживающих методов, по зашумленным данным оцениваются функция и ее производные, а затем на втором этапе строятся МНК-оценки параметров уравнения. Двухэтапный метод легко реализуется, однако он может быть статистически неэффективным, так как производные не могут быть точно оценены по зашумленным данным, особенно производные старших порядков. Предложенная методика верификации основана на сопоставительном анализе параметров Рисунок 1 Временной ряд и результаты авторегрессионной модели и разностного уравнения. моделирования Полученные результаты показывают возможность Экспериментальным подтверждением использования такого подхода при моделировании правильности предложенной методики должно быть зашумленных динамических объектов с хорошее совпадение наблюдаемых значений распределенными переменными и обосновывают температуры, значений, полученных по проведение дальнейших исследований в этом авторегрессионной модели (10), и численного направлении. решения уравнения диффузии. Качество авторегрессионной модели и модели диффузии Литература характеризуется коэффициентом детерминации 𝑅2 . [1] Захаров, Е.В., Дмитриева И.В., Орлик С.И.: В обоих случаях этот коэффициент равен 0,78, что Уравнения математической физики. говорит об одинаковых качественных Университетский учебник. М.: Издательский характеристиках моделей. Степень близости центр «Академия», 320 с. (2010) решений, полученных по этим моделям можно [2] Dickey, D.A., Fuller, W.A.: Distribution of the оценить коэффициентом парной корреляции, Estimators for Autoregressive Time Series with a который составил в нашем случае 0,97, что еще раз Unit Root. J. of the American Statistical подтверждает практическую эквивалентность Association, 74 (366), pp. 427-431 (1979). doi: найденных решений. 10.2307/2286348 5 Заключение [3] Носко, В.Н.: Эконометрика. Кн. 2. Ч. 3, 4. Учебник. М.: Изд. дом «Дело» РАНХиГС, 576 с. Для идентификации динамических моделей с (2011) распределенными переменными обычно [4] Мареев, В.В., Станкова, Е.Н.: Основы методов используются либо непараметрические методы конечных разностей. СПб.: Изд-во С.-Петерб. идентификации, например, на основе результатов ун-та, 64 с. (2012) активного эксперимента [11], который не всегда возможен, либо параметрические методы с [5] Clements, M.P., Hendry, D.F.: Forecasting Non- аппроксимацией производных [12]. В последнем stationary Economic Time series. Cambridge, случае возникает опасность неадекватного Massachusetts, London: MIT Press, 262 p. (1999) 357 [6] Patterson, K.: An Introduction to Applied American Statistical Association. 108 (503), Econometrics: A Time Series Approach. New pp. 1009-1020 (2013). doi: York: Palgrave, 832 p. (2000) 10.1080/01621459.2013.794730 [7] Канторович, Г.Г.: Анализ временных рядов. [14] Putter, H., Lange S.: A Bayesian Approach to Экономический журнал ВШЭ. (3), сс. 379-401 Parameter Estimation in HIV Dynamical Models. (2003) Statistics in Medicine Heisterkamp, pp. 2199-2214 [8] Хромов, С.П., Петросянц, М.А.: Метеорология и (2000) Климатология. Учебник. М.: Изд-во МГУ, 527 с. [15] Huang, Y., Liu, D.: Hierarchical Bayesian Methods (2001) for Estimation of Parameters in a Longitudinal HIV [9] Матвеев, М.Г., Михайлов, В.В., Сирота Е.А.: Dynamic System, pp. 413-423 (2006) Комбинированная прогностическая модель [16] Ramsay, J. O., Hooker, G., Campbell, D.: нестационарного многомерного временного Parameter Estimation for Differential Equations: a ряда для построения пространственного Generalized Smoothing Approach (with профиля атмосферной температуры. discussion). J. of the Royal Statistical Society, Информационные технологии, 22 (2), сс. 89-94 Series B, pp. 741-796 (2007) (2016) [17] Cao, J., Huang, J., Wu, H.: Penalized Nonlinear [10] NCEP/DOEAMIPIIReanalysis [Электронный Least Squares Estimation of Time-varying ресурс]. http://www.cdc.noaa.gov/cdc/data.ncep. Parameters in Ordinary Differential Equations. J. of reanalysis2.html Computational and Graphical Statistics, pp. 42-56 [11] Бойков, И.В., Кривулин, Н.П.: Методы (2012) идентификации динамических систем. [18] Liang, H., Wu, H.: Parameter Estimation for Программные системы: теория и приложения. Differential Equation Models Using a Framework 23 (5), С. 79-96 (2011) of Measurement Error in Regression Models. J. of [12] Bar, M., Hegger, R, Kantz, H.: Fitting Differential the American Statistical Association, pp. 1570- Equations to Space-time Dynamics. Physical 1583 (2008) Review, E. 59 (1), pp. 337-342 (1999). doi: [19] Chen, J., Wu, H.: Efficient Local Estimation for 10.1103/PhysRevE.59.337 Time-varying Coefficients in Deterministic [13] Xun, X., Cao, J., Maity, J.: Parameter Estimation of Dynamic Models with Applications to HIV-1 Partial Differential Equation Models. J. of the Dynamics. J. of the American Statistical Association, pp. 369-384 (2008) 358