УДК 681.2.088+519.254 Хрусталёв Юрий Петрович,
к. т. н., доцент кафедры вычислительной техники ИрГТУ, тел. 40-51-07, e-mail: [email protected]
Акулов Владислав Михайлович,
начальник отдела ВСФ ВНИИФТРИ, учёный хранитель Вторичного эталона времени, частоты и шкалы времени,
тел. 46-73-30, e-mail: [email protected] Ипполитов Александр Александрович, аспирант кафедры вычислительной техники ИрГТУ, тел.: 40-51-07, e-mail: [email protected]
Курышева Людмила Николаевна,
начальник метрологического сектора ВСФ ВНИИФТРИ, тел. 46-73-30, e-mail: [email protected]
ПОВЫШЕНИЕ ТОЧНОСТИ ОЦЕНИВАНИЯ ОТНОСИТЕЛЬНЫХ ОТКЛОНЕНИЙ ЧАСТОТЫ ГРУППОВОГО ЭТАЛОНА
Y.P. Khrustalev, V.M. Akulov, A.A. Ippolitov, L.N. Kurysheva
INCREASING ACCURACY OF RELATIVE FREQUENCY DEVIATIONS ESTIMATION IN COLLECTIVE STANDARD
Аннотация. Для обеспечения единства измерения единиц времени и частоты на территории Российской Федерации действует Государственная служба времени, частоты и определения параметров вращения Земли (ГСВЧ РФ), основу которой составляют Государственный эталон и ряд вторичных эталонов. Вторичный эталон ВЭТ1-5 функционирует в составе ВСФ ВНИИФТРИ (г. Иркутск). В работе рассматриваются способы повышения стабильности эталона ВЭТ1-5, связанные с совершенствованием методов обработки измерительной информации.
Ключевые слова: групповые эталоны, стабильность, статистика, модели, АРПСС.
Abstract. State Service for Time and Frequency (SSTF) of Russia operates with the aim of providing uniformity of time and frequency measurements across the whole territory of country. It is based on National standard and a number of secondary standards. The secondary standard VET1-5 operates in East-Siberian Branch of VNIIFTRI (Irkutsk). Ways of increasing stability of Secondary standard for time, frequency and time scale VET1-5 related to enhancement of measuring data processing methods are considered in this work.
Keywords: collective standards, stability, statistics, models, ARIMA.
В настоящее время наиболее точно воспроизводимой физической величиной являются единицы времени и частоты. Хранение этих величин обеспечивается Государственным эталоном и вторичными эталонами. Для доведения точностных характеристик Государственного эталона до потребителей создана сеть вторичных эталонов и
система передачи размеров единиц времени и частоты и шкалы времени первичного эталона иТС (8И) рабочим эталонам и рабочим средствам измерений. Ответственность за обеспечение единства время-частотных измерений на территории страны возложена на Государственную службу времени, частоты и определения параметров вращения Земли (ГСВЧ РФ).
Вторичный эталон ВЭТ1-5, функционирующий в составе Восточно-сибирского филиала ВНИИФТРИ, является средством высокоточного хранения размеров единиц времени и частоты, а также шкалы времени иТС ^И) первичного эталона с суммарной погрешностью, не превышающей ± 2 • 10 14 по частоте и ± 80 наносекунд по шкале времени. Вторичный эталон ВЭТ1 -5 является одним из элементов группового эталона ГСВЧ, хранители времени и частоты которого используются при формировании групповой шкалы системы ГСВЧ.
Ведение вторичного эталона в системе ГСВЧ с указанными выше метрологическими характеристиками осуществляется на базе сложных технических и структурно взаимосвязанных аппаратно-программных комплексов.
Основным аппаратурно-программным комплексом вторичного эталона ВЭТ1-5 является аппаратура хранения размеров единиц частоты и времени. В настоящее время эта аппаратура состоит из шести водородных стандартов частоты и времени типа Ч1-75А (четыре прибора), Ч1-70М (один прибор), Ч1-70 (один прибор) с системами автоматической настройки резонаторов. Лучшие приборы этого комплекса характеризуются суточной нестабильностью частоты менее 1 • 10"15.
ш
Не менее важными аппаратно-программными комплексами, образующими в целом структуру, именуемую вторичным эталоном времени и частоты, являются:
• система внутренних сличений эталона,
• система внешних сличений эталона,
• система формирования сигналов физических шкал времени,
• управляющие программные средства систем внутренних и внешних сличений,
• алгоритмы и программные средства, используемые для формирования аналитических и физических шкал времени шкал времени эталона,
• система буферизации и размножения сигналов времени и частоты,
• аппаратура и программные средства системы информационного обмена в рамках локальной вычислительной сети, а также внешний БТР-сервер,
• аппаратурный комплекс системы жизнеобеспечения эталона ВЭТ1-5.
В настоящей работе рассматриваются вопросы обработки данных, полученных в системе взаимных измерений эталона. Проблема повышения стабильности групповых эталонов за счёт уменьшения методической погрешности оценивания относительных отклонений частоты опорного элемента группового эталона весьма актуальна, поскольку эти результаты применимы и для приборов другого типа. Т. е. при повышении метрологических характеристик стандартов частоты и времени, которое может быть достигнуто за счёт разработки новых типов стандартов (что является весьма дорогостоящим мероприятием), точность группового эталона может быть дополнительно увеличена ещё на несколько процентов за счёт применения соответствующих алгоритмов [1].
При формировании автономной шкалы времени группового эталона производится измерение разностей относительных отклонений частот А///0 стандартов, входящих в состав эталона. При этом один из стандартов выбирается в качестве опорного (ведущего) и измеряется разность
АГ"//0оп-А/7/, . = 1,...,п,
где п - число стандартов, входящих в эталон.
Обозначив А/ //о как у, а А/оп //°п -как уоп, получаем систему линейных уравнений
Ъ = Уоп -У. * =1 п (1)
Для определённости будем считать опорным стандарт с номером 1, что не меняет общности рассуждений. Чтобы упростить запись, дополним систему (1) «фиктивным измерением»
= у - Л =0 .
Мы имеем недоопределённую систему уравнений, в которой число уравнений на 1 меньше числа неизвестных. (Погрешностью измерений можно пренебречь, что вполне допустимо при обработке данных на суточных интервалах. ) Таким образом, задача формирования автономной шкалы времени сводится к решению недоопределённой системы уравнений (1), поскольку при известных относительных отклонениях частоты стандартов А///0 поправка к показаниям часов на интервале времени т находится из соотношения А///о ' т.
Известно, что недоопределённые системы (или, в случае линейных систем, системы с неполной матрицей наблюдений) не имеют единственного решения. Чтобы получить такое решение, необходимо ввести в систему дополнительные ограничения. Чаще всего в качестве ограничения принимают условие минимальной нормы вектора решения [2]. В этом случае решение может быть найдено через псевдообратную матрицу. При существующей системе измерений это решение даёт значение относительного отклонения частоты опорного элемента (стандарта частоты) как [3]:
1 п
•Рои = - 'Е (2)
п .=-
(в выражении (2) считаем, что = уоп -уоп = 0).
Значения относительных отклонений частоты остальных элементов легко находятся из соотношений (1) и (2) либо непосредственно через псевдообратную матрицу.
Нетрудно показать, что решение (2) можно было бы получить, введя в систему уравнений (1) вместо фиктивного соотношения =0 условие
п
ЕУ = 0 . Анализ погрешности решения [3] пока-
.=1
зывает, что значение частоты опорного стандарта, найденное по формуле (2), соответствует «истинному» её значению только при выполнении последнего условия, т. е. сумма относительных отклонений частоты элементов группового эталона должна равняться нулю. Разумеется, это весьма жёсткое условие и вряд ли оно может быть выполнено при небольшом числе стандартов, входящих в групповой эталон.
Точность оценивания относительного отклонения опорного элемента группового эталона (стандарта частоты) можно увеличить, если помимо результатов измерений на каждом такте измерений (в дальнейшем будем иметь в виду измерения на суточных интервалах) учитывать ещё и прогнозы этих отклонений для каждого из элемен-
тов группового эталона. Формула, используемая для оценки Aj/J0on , получена при декомпозиции фильтра Калмана [4] и имеет следующий вид:
Д^ оп n n
J = Уоп = Е [Уоп + (У - У )] = Е [Z + y (1)] ' (3) Jo 1=1 1=1
где z = Уоп - У - результат измерения относительных разностей частот опорного и i-го стандарта, g - вес измерения,
У (1) - прогноз относительного отклонения частоты /-го генератора, вычисленный на предыдущем такте обработки данных (с упреждением на один шаг).
Погрешность оценивания по формуле (3)
равна
А = Ё& [У - У (1)]
(4)
Очевидно преимущество использования формулы (3) по сравнению с выражением (2), поскольку теперь равенство нулю погрешностей оценивания может быть достигнуто и без выпол-
п
нения условия ^у = 0, т. е. условие (3) - более
i=1
«мягкое» условие, заключающееся в равенстве нулю взвешенной суммы отклонений прогнозов у от их истинных значений. (Напомним, что речь идёт об обработке данных на текущем такте 1. Веса g¡ выбираются обратно пропорциональными дисперсии прогнозов.)
Полученные результаты согласуются с предложением Ж. Рютмана [5] рассматривать нестабильность как степень непредсказуемости значений частоты. Действительно, как бы ни изменялась частота высокостабильных генераторов, если мы знаем закон её изменения, всегда можно ввести соответствующую поправку.
Таким образом, задача оценивания относительных отклонений частоты опорного элемента группового эталона сводится, помимо измерений взаимных разностей частот, к получению краткосрочных прогнозов у (1) .
Для вычисления прогнозов у (1) необходимо иметь математические модели рассматриваемых процессов. Поскольку речь идёт о значениях частот (относительных отклонений частот) для равноотстоящих моментов времени, целесообразно использовать соответствующие модели - модели временных рядов. В качестве таких моделей широко используются динамические стохастические модели - модели авторегрессии - проинтегрированного скользящего среднего (АРПСС).
Преимущество использования таких моделей перед полиномиальными моделями показано в [6]. В настоящее время разработана методика построения моделей АРПСС, имеются пакеты прикладных программ, реализующие данную методику. Однако для её использования необходимо иметь исходные временные ряды. Как указывалось выше, результаты измерений, выполняемых в процессе ведения эталона времени и частоты, представляют собой разности частот опорного и /-го стандарта. Т. е. мы имеем дело с результатами косвенных измерений и напрямую применить методику Бокса -Дженкинса для построения моделей АРПСС невозможно. В [3] предложено для решения поставленной задачи минимизировать сумму квадратов результатов измерений, выполняемых на каждом такте обработки данных, и их прогнозов. Поскольку прогнозы измерений зависят от вектора параметров используемых моделей, то оптимальные оценки этого вектора соответствуют минимуму функционала. При известной структуре моделей (т. е. при заданном числе членов авторегрессиир и скользящего среднего д) процесс нахождения оптимальных параметров сходится. Если информация о структуре модели отсутствует, можно воспользоваться методом простого перебора, последовательно наращивая величины р и д для каждого из стандартов. Т. к. числа р и д обычно невелики (не превышают 5 для большинства реальных систем [6]), то число рассматриваемых структур конечно и можно рассмотреть все возможные варианты и выбрать модель с минимальным значением функционала.
В [1] рассмотрен другой подход к решению данной проблемы. Исходя из заданной матрицы измерений, находим значения относительных отклонений частоты опорного стандарта для каждого момента времени X через псевдообратную матрицу. Результат находится как среднее арифметическое всех измерений. Затем вычисляем соответствующие характеристики всех остальных элементов группового эталона. Получаем временные ряды - оценки относительных отклонений частоты , для которых находим порядки авторегрессии р^ и скользящего среднего дг-. Применяя методику Бокса - Дженкинса, строим соответствующие модели АРПСС, параметры которых используем в качестве начальных оценок в процедуре оптимизации, описанной выше.
Предложенная методика позволила уменьшить погрешность воспроизведения частоты вторичного эталона ВЭТ1-5 приблизительно на 8 % (по результатам внешних сличений с первичным эталоном на трёхмесячном интервале). Кажущееся незначительным снижение погрешности воспро-
1=1
ш
изведения частоты вторичного эталона на самом деле может привести к существенному экономическому эффекту при учёте высокой стоимости аппаратуры.
Особого внимания заслуживает проблема учёта детерминированных трендов частоты, имеющих место в реальных системах. Тренды могут быть учтены различными способами:
• введением члена 0О в модели авторегрессии проинтегрированного скользящего среднего, построенные для разностных рядов;
• непосредственным введением детерминированных функций в исходные временные ряды.
Против первого подхода резко возражают такие авторитетные источники, как [7]. Нам также представляется более предпочтительным второй подход, т. к. в этом случае параметры используемых детерминированных функций имеют ясное физическое толкование. Собственно, анализ результатов измерений при статистической обработке данных, а сейчас мы говорим именно об обработке данных в режиме их накопления, должен начинаться с обнаружения детерминированных функций во временных рядах и оценивания их параметров.
Поскольку мы имеем дело с недоопределён-ной системой, однозначно оценить параметры детерминированных трендов по результатам измерений, выполняемых в эталоне, невозможно. Чтобы сохранить «привязку» вторичного эталона к первичному, следует использовать результаты внешних сличений, из которых можно найти оценки параметров тренда опорного стандарта частоты. Определение параметров трендов остальных элементов по результатам взаимных измерений после этого не представляет принципиальных трудностей. Чаще всего имеют место линейные тренды частоты, хотя возможны и тренды более высоких порядков.
Если по какой-либо причине мы не имеем возможности (или не хотим ею воспользоваться) привлечь при анализе данных результаты внешних сличений эталона, необходимо проверить гипотезу об отсутствии линейной функции в рядах измерений ъ{(), ■ = 2,.,п . Это можно достаточно просто сделать, используя специализированные пакеты прикладных программ, предназначенные для статистического анализа. Если хотя бы один ряд )
не содержит линейных функций, можно с большой долей вероятности считать, что линейный тренд отсутствует в рядах, соответствующих опорному и у-му стандарту. Оценивание параметров линейных функций в остальных рядах - задача тривиальная. Если же гипотеза об отсутствии линейного тренда хотя бы в одном ряде не подтверждается, придётся
полагать, что он отсутствует у временного ряда с минимальным углом наклона линейной функции.
После того, как произведена оценка линейных составляющих в рядах относительных отклонений частоты у (), следует устранить их влияние и произвести построение моделей авторегрессии скользящего среднего (АРСС) в соответствии с описанной выше методикой.
Изложенные процедуры представляют собой статистическую обработку результатов измерений, выполняемых в процессе ведения группового эталона. В результате получаются ряды относительных отклонений частот у, (т. е. А/1 //. ) для каждого из стандартов, входящих в состав эталона, а также математические модели этих рядов: оценки углов наклона линейных трендов и параметры моделей АРСС, т. е. коэффициенты фг, 9г и оценка
остаточной дисперсии о а. Наличие математических моделей позволяет решать задачи динамической обработки данных, т. е. задачи фильтрации.
Под динамической обработкой будем понимать задачи получения оценок вектора состояния
эталона Ут =(у,..., у) по результатам измерения
2т = (ъ ,.■■, ) в темпе их поступления. При этом равенство размерности векторов У и Z достигается за счёт введения фиктивного измерения = 0, как это делалось выше. Задачи статической и динамической обработки детально рассматривались в [3]. При этом динамическая обработка основывается на использовании соотношения (3), где прогнозы вычисляются на основе ранее построенных математических моделей. Прогнозы содержат детерминированную составляющую (учёт тренда) и составляющую, вычисленную с использованием моделей АРСС.
В настоящей работе рассмотрим особенности динамической обработки данных, обусловленные характеристиками водородных стандартов частоты. Дело в том, что ряды относительных отклонений частоты водородных стандартов помимо детерминированных трендов и составляющих, описываемых уравнениями АРСС, могут содержать ещё и ступенчатые функции, представляющие собой скачки частоты. Моменты и амплитуды скачков заранее не известны, и их идентификация и учёт непременно должны выполняться в процессе динамической обработки данных.
Оценки о2а1, а следовательно, и оа1, полученные в процессе построения моделей АРСС, позволяют достаточно просто решить эту задачу. С помощью построенных моделей вычисляются прогнозы у (1) для каждого из рядов оценок относительных отклонений частоты от их доверительной границы. В частности, при доверительных
уровнях 0,99 прогнозы лежат в границах ± 3оа. После вычисления оценок у (^) на текущем такте
обработки данных I, проверяются условия попадания оценок в доверительные интервалы. Если все оценки попадают в эти интервалы, считается, что на данном такте скачков частоты не обнаружено. Если для какого-либо из стандартов частоты это условие не выполняется, данный стандарт исключается на текущем шаге из обработки и вновь вычисляются оценки у (^) (но уже без учёта результата измерения, содержащего стандарт, подверженный скачку частоты). Оценка у (^) находится
из результата измерения (^) и полученной оценки у (^) - относительного отклонения частоты опорного стандарта. Если после исключения из обработки измерений, содержащих «скачок частоты», и вычисления оценок у{) для оставшихся вновь идентифицируется скачок, процедура повторяется и т. д. до тех пор, пока все оценки уг (^)
не будут находиться в доверительных интервалах. Т. е. строится итерационная процедура отбраковки «сомнительных» измерений. Если на каком-либо шаге итерационной процедуры выходят за доверительные границы несколько значений yj ), из
обработки исключается измерение, соответствующее стандарту, частота которого наиболее далеко вышла за доверительные границы. Однако при этом рассматривают не просто отклонение частоты от границ прогнозов, а их значения, нормированные оценками с>г , т. е. рассматриваются вели-
чины
Ау.
и отбрасывается результат измерения
Рис. 1. Упрощённый алгоритм отбраковки скачков частоты
Существует много критериев контроля адекватности [6]. Для наших целей представляется наиболее удобным использовать статистику:
Е <
п+к
(5)
Е аа + Е <
г! (^), соответствующий максимуму этих величин.
При выходе всех оценок за доверительные границы считается, что имел место скачок частоты опорного стандарта.
Упрощённый алгоритм отбраковки скачков частоты имеет вид, приведенный на рис. 1.
Весьма важным при вычислении прогнозов является контроль адекватности прогнозирующих моделей. Речь идёт об обнаружении моментов «разладки» моделей. На самом деле происходит не разладка моделей, а изменение статистических характеристик временных рядов, описывающих относительные отклонения частоты стандартов. Так или иначе, модели, используемые для прогнозирования, становятся неадекватными реальным процессам. Необходимо произвести их перенастройку.
] = п+1
где а - остатки, вычисленные в процессе построения модели,
п - объём выборки, по которой подгоняются параметры модели,
а ■ - остатки от прогнозов, используемых в процессе функционирования алгоритма, к - текущий шаг обработки данных.
^-статистика позволяет сопоставлять сумму квадратов остаточных членов, полученную на этапе построения модели, с текущим её значением и имеет функцию плотности вероятности [8]:
п + к Л 2 "
Г|
/ (х) =
Г
пЧ к
-• х
1.(1 - X )а-1
(6)
2 ) \ где Г - гамма-функция.
2
1=1
п
2
1=1
ш
Функция /(х), определяемая выражением (6), является частным случаем бета-распределения:
•(1 - х)■"', (7)
п к
при р = п2; я=2
При ухудшении качества прогнозов, т. е. при «уходе» процесса от ранее построенной модели, вклад второй суммы знаменателя в выражении (5) начинает возрастать, что повлечёт уменьшение величины X. Если при заданном уровне значимости величина X, вычисленная по формуле (5), ока-
п
жется меньше значения р с параметрами р = — ; к
Я = —, то математическую модель следует считать
неадекватной реальному процессу, т. е. модель нуждается в перенастройке.
Разумеется, при использовании данного критерия, как, впрочем, и остальных, основанных на предположении нормальности остатков, необходимо провести проверку его выполнения.
Для повышения эффективности алгоритма динамической обработки данных следовало бы рассмотреть процедуру адаптации моделей, производя их подстройку пропорционально остаткам от прогнозов, однако это выходит за рамки настоящей работы.
Выводы
1. Точность воспроизведения частоты групповыми эталонами времени и частоты можно повысить, используя алгоритмы оценивания относительных отклонений частоты, основанные на использовании их прогнозов. При этом выполняется предположение Ж. Рютмана об оценке нестабильности частоты как меры её непредсказуемости.
2. В процессе статической обработки данных (обработки данных в процессе их накопления) возможно построение математических моделей рядов относительных отклонений частоты по результатам косвенных измерений, т. е. при отсутствии исходных временных рядов. Процедура построения моделей АРСС основана на минимиза-
ции суммы квадратов отклонений результатов измерений от их прогнозов.
3. При функционировании эталонов целесообразно применять процедуру динамической обработки данных (процедуру фильтрации), решая при этом проблему отбраковки аномальных измерений (выбросов).
4. Предлагаемые в работе алгоритмы могут быть использованы при построении групповой шкалы времени территориально разнесённых локальных эталонов времени и частоты. При этом необходимо учитывать шумы измерений, выполняемых с помощью системы ОР8/ГЛОНАСС.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Хрусталёв Ю. П., Акулов В. М., Ипполитов А. А., Курышева Л. Н. Обработка данных, полученных по результатам взаимных измерений вторичного эталона времени и частоты // Вестник ИрГТУ. 2012. № 7. С. 22-28.
2. Гантмахер Ф. Р. Теория матриц. М. : Наука, 1966. 576 с.
3. Хрусталёв Ю. П. Статическая и динамическая обработка данных, получаемых в процессе ведения эталонов времени частоты // Измерительная техника. 2004. № 6. С. 20.
4. Ипполитов А. А., Хрусталёв Ю. П. Субоптимальная фильтрация в системах с неполной матрицей наблюдений // Информационные и математические технологии в науке и управлении : тр. XV Байкал. Всерос. конф. Ч. 1. Иркутск : ИСЭМ СО РАН, 2010. С. 174-182.
5. Рютман Ж. Характеристики нестабильности фазы и частоты сигналов высокостабильных генераторов. Итоги развития за 15 лет // ТИИ-ЭР. 1978. Т. 66. № 9. С. 70-102.
6. Бокс Д., Дженкинс Г. Анализ временных рядов. Прогноз и управление. Вып. 1. М. : Мир, 1974. 406 с.
7. Тюрин Ю. Н., Макаров А. А. Статистический анализ данных на компьютере. М. : Инфра-М, 2003. 544с.
8. Крамер Г. Математические методы статистики. М. : Мир, 1975. 648 с.