Фрактальный анализ состояния неустойчивости физической системы Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

УДК 004.42:539.1

Паньгин Александр Викторович

ФРАКТАЛЬНЫЙ АНАЛИЗ СОСТОЯНИЯ НЕУСТОЙЧИВОСТИ ФИЗИЧЕСКОЙ СИСТЕМЫ

Аннотация

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

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

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

Джефф О'Лири

ВВЕДЕНИЕ

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

© Паньгин A.B., 2012

ния сигнала могут находиться на границе допустимой неопределенности (например, для средневзвешенной мощности реактора ВВЭР данный показатель составляет ~ 2 %, [1]).

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

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

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

Для интерпретации регистрируемых сигналов нейтронно-физических датчиков с целью оперативного контроля и обеспечения условий безопасной эксплуатации РУ разработаны эффективные методы виброшумовой диагностики [2].

В работе [3] показано, что наличие флуктуаций нейтронного потока в характерном диапазоне частот 0.1-3.0 Гц, вызванных кипением, может быть зарегистрировано датчиками системы внутриреакторного контроля (СВРК).

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

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

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

- изменением характеристик измери-

тельных трактов (первичных преобразователей);

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

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

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

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

Для анализа сложных динамических систем все большее распространение находят представления фрактальной геометрии [4], отображающие структуру исследуемых процессов реальных нелинейных самоорганизующихся систем. Критерий фрактальности сигнала может служить дополнительным диагностическим фактором, характеризующим состояние физической системы локально и в целом. Методика определения фрактальных характеристик размерности временных рядов используется во многих областях при решении задач как идентификации (определения типа поведения), так и для прогнозирования различных процессов [5].

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

параметров) с помощью локального индекса вариации для диагностирования стохастической природы исследуемого сигнала в режиме теплогидравлической неустойчивости.

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

МЕТОДИКА ИССЛЕДОВАНИЯ

Основной характеристикой компактного фрактального множества А является его размерность Б/, определяемая в прикладных задачах для произвольного евклидового метрического пространства следующим образом:

1 5®0 1п(1/ 5)

где N^5) - минимальное число примитивов (например, шаров) с геометрическим фактором (линейным размером) 5, образующих покрытие множества А. Показано [4], что в эквивалентных метриках размерность множества А, определяемая формулой (1), не меняется, а следовательно, для ее определения может быть использовано клеточное покрытие и линейные преобразования координатных масштабов.

Формулу (1) можно заменить эквивалентной

^(5) ~ (1/5)°/. (2)

Поскольку для покрытия единичного отрезка, квадрата или куба его копиями размера 5, их потребуется соответственно 1/5, 1/52 и 1/53, то показатель Б/ является аналогом топологической размерности °т, но в отличие от последнего может быть дробным числом, при этом величина / = - °т определяется как индекс фрактальности.

Вычисление по формуле (1) практически затруднено, так как слишком медленно

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

Фрактальность временного ряда может быть охарактеризована с помощью показателя Херста (Hurst), который вычисляется, например [6], следующим образом.

Пусть x (t) - анализируемый участок временного ряда длины N. Тогда для различных вложенных (целиком в длину N) периодов t определяется взаимосвязь

(R/S)^°< R/S> ~ctH, (3) где с, H - некоторые константы; <...> - усреднение по числу (m = N/t) укладывающихся периодов длины t во временной участок ряда длины N; Rt , St - величины, рассчитываемые для каждого периода t по формулам:

t t Rt = max У (x(u) - Xt) - min У (x(u) - xt) -

1£t£t '

u=1 u=1

размах накопленных сумм временного ряда;

= 1У (х(и) - х. )2 - среднеквадратичное отклонение ряда; X. =1 ^ х(и) - сред. и =1

нее значение ряда.

Константа Н (показатель Херста) характеризует угол наклона графика линейной регрессии зависимости 1п(Л/£)т от 1и(т) и связана с формулой: Б/= 2 - Н. Значение показателя Херста может изменяться от 0 до 1 и позволяет определить при подходящем масштабировании нормированную меру корреляции С(?) временных событий х (?) и х (-?), полагая, что х(0) = 0, [6]:

С(,) = < х(-)) > = 22н-1 -1. (4)

< х(0 >

При Н> 0.5 корреляция положительна, то есть тенденция к увеличению в прошлом означает в среднем тенденцию к увеличению

в будущем, и наоборот, тенденция к уменьшению в прошлом означает тенденцию к уменьшению в будущем. Такой процесс называется персистентным (сохраняющим тенденцию).

При Н < 0.5 корреляция отрицательна, то есть увеличение в прошлом означает вероятное уменьшение в будущем, а тенденция к уменьшению в прошлом означает увеличение в будущем. Такой процесс характеризуются отсутствием устойчивости (анти-персистентностью).

Значения Н, близкие к 0.5, характеризуют неопределенность поведения (некоррелированность отсчетов) типа винеровского случайного процесса (модель классического броуновского движения, см. Приложение).

Описанный статистический метод называется также методом нормированного размаха или Л/^-анализом. Однако для надежного вычисления значения Н требуется большая репрезентативная выборка, например, алгоритмом программного кода БЯЛСТАМ (http://www.impb.ru/downloads/ fractan.zip) производится анализ временного ряда, содержащего несколько тысяч данных. Внутри такого масштаба временной ряд может менять характер своего поведения (например, проявляются различного рода флуктуации сигнала).

Чтобы связать локальную динамику соответствующего процесса с фрактальной размерностью ряда на некотором временном интервале, на практике определяется «промежуточная асимптотика» размерности (то есть на масштабах, меньших характерного времени изменения основных динамических состояний процесса, но больших минимального масштаба 80, связанного с периодом опроса измеряемого параметра).

В работе [5] локальная фрактальная структура временного ряда исследуется с помощью индекса вариации, который определяется по изменению амплитуды сигнала в «скользящем окне» следующим образом.

Выбирается «скользящее окно» (О из п +1 отсчетов временной длительностью [г0, гп], которое равномерно разбивается на т отрезков размером 8=(гп- г0)/т. На каж-

дом отрезке [^ ^ г8], где I = 1,2, ...,т, определяется размах А1 (8) отсчетов (разность между максимальным и минимальным значениями ряда х(г) на этом временном отрезке).

Рассчитывается значение вариации сигнала на временном интервале (

т

Уа (8) = £ А, (8), (5)

,=1

которое имеет (для различных 8) подходящее представление:

Уа(8) ~ 8-". (6)

Показатель V характеризует угол наклона линейного тренда графика точек (1п (8), 1п( У()) - значений логарифмов длин интервалов и логарифмов вариации для различных разбиений.

Индекс V, как следует из формулы (2) связан в пределе с фрактальной размерностью формулой: V = и совпадает с индексом фрактальности ¡л для одномерных структур (Бт = 1). Соответственно интерпретация индекса V, как показателя устойчивости, согласуется с представлениями для показателя Херста из соотношения: Н = 1- V [5].

Задача идентификации макросостояния физической системы может быть сведена к отслеживанию значения «фрактального» индикатора наблюдаемого временного ряда показаний следящей измерительной аппаратуры.

РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЯ

Апробирование методики проведено на ряде характерных экспериментов, зарегистрированных ИВК в режимах с заведомым переходом через границу устойчивости.

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

тракта ИВК, амплитуда и период которых существенно меньше, чем у колебаний сигнала в режиме неустойчивости физической системы.

Для выделения полезного сигнала и подавления неинформативного шума производилась предварительная обработка данных с помощью вейвлет-преобразования средствами Wavelet Toolbox пакета Matlab. Для сигналов с локальными особенностями вейвлет-преобразование более эффективно, чем использование рядов Фурье [7].

Суть вейвлет-преобразования одномерного сигнала x (t) состояла в его разложении (7) по базисной функции (8) - вейвлету Y(t) посредством ее масштабных изменений и сдвигов (9).

-¡F k (t )■

(7)

x(t ) = X C

j ,k=-¥

Выбирался простейший HAAR-вейвлет, определяемый выражением:

1, 0 < t < 1/2; ¥ (t ) = <-1, 1/2 < t < 1; (8)

0, t < 0, t > 1.

Ортогональные базисные функции имеют вид:

Yjk(t) = 2 J/2Y(2 jt - k). (9)

Для выделения «полезного» сигнала использовались превалирующие коэффициенты вейвлет-разложения исходного сигнала нейтронной мощности с помощью наглядной оболочки-меню «Wavelet Coefficients Selection 1-D».

Последующий фрактальный анализ «синтезированного» временного ряда проводился средствами MS Excel. При вычислении индекса вариации n достаточно использовать, как показано в [5], последовательность небольшого числа вложенных разбиений окна w. Для практической реализации выбирались m = 2k, где k = 0, 1, ..., 5. Чтобы соотнести значение n (t) с поведением временного ряда в точке t, определялось значение n на интервале, предшествующем точке t (при анализе дискретного отклика датчика вы-

биралось окно, состоящее из 33 временных отсчетов).

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

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

На рис. 1 представлен график показаний сигнала мощности в эксперименте 1, «очищенного» от шума, в «скользящем окне» из 33 временных отсчетов (предшествующих некоторой временной точке г) в режиме теп-логидравлической неустойчивости течения теплоносителя.

Отметим, что флуктуации сигнала могут иметь периодический или нерегулярный характер. Из графика видно, что при неустойчивом режиме, колебание мощности активной зоны может достигать более 3% N

^ ном

(номинальной мощности РУ).

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

Рис. 1. Фрагмент амплитудных колебаний корректированного (кор.) сигнала мощности в режиме теплогидравлической неустойчивости

3,5 3,0 2,5

г,о

1.5 1,0 0,5 0,0

<к---«-..

ч = -О.бОВу + 3 К2 = 0.998 288

12 3 4 ■Вариация сигнала;-Тренд

Рис. 2. Результат вычисления индекса вариации

щью метода наименьших квадратов (МНК) линию регрессии (тренд) у = ах + Ь и отождествить V = -а.

На рис. 3 представлены результаты обработки для эксперимента 1 временного ряда сигнала нейтронной мощности % N (в процентах от номинального значения, регистрируемого измерительным комплексом при испытании активной зоны реактора) для определения показателя V (индекса вариации) его фрактальной структуры.

Результаты исследования методом фрактального анализа временного ряда сигнала мощности, зарегистрированном в данном эксперименте, показали, что индекс вариации V (?) ряда изменяется от значения V ~ 0, что соответствует топологической размерности (Су = Ст = 1) графика временной зависимости значения сигнала при «нормаль-

Рис. 3. Идентификация режима неустойчивости сигнала мощности РУ методом фрактального анализа

ном» режиме работы ЯЭУ, до значения V ~ 0.4-0.7, что соответствует фрактальной размерности (Су > 1.4) сигнала и хаотичности характера изменения его амплитуды в режиме теплогидравлической неустойчивости, а также позволяет определить временные отсчеты границы режима неустойчивости в работе установки от начала эксперимента.

Примечание: фрактальный индикатор в переходных устойчивых маневрах мощности (подъема, спада) мал по величине, однако чувствителен к локализованным одиночным выбросам. На рис. 3 одиночный скачок сигнала мощности (А N = 1.2 %) вызвал изменение V~ 0.3.

На рис. 4 приведен фрагмент графика изменения значения индекса вариации (V) временного ряда, соответствующего регистрации нейтронной мощности (% N ) в эксперименте 2 после предварительной обработки удаления «шумовой» компоненты исходного сигнала при помощи вейвлет-ана-лизатора.

Как и для ранее рассмотренного эксперимента 1 проявление признаков теплогид-равлической неустойчивости состояния активной зоны характеризуется соответствующим значением индекса вариации V > 0.4. Однако для данного эксперимента режим неустойчивости имеет более сложный характер с течением времени. Локальная фрактальная размерность ряда изменяется от 1 до 1.87. Развернутый анализ теп-логидравлической неустойчивости выходит за рамки формата исследований и потому не приводится.

ЗАКЛЮЧЕНИЕ

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

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

Показана возможность эффективного применения фрактального анализа сигналов в задаче идентификации неустойчивых режимов работы реакторной установки. На основе численных данных алгоритмов обработки установлена эмпирическая зависимость между локальным значением индекса вариации на выбранном «скользящем» временном интервале и детерминированностью процесса на этом периоде. Отмечено, что величина индекса вариации V > 0.4 определяет появление флуктуа-ций сигнала с признаками неустойчивости. Используемые алгоритмы фрактального анализа достаточно просты и эффективны в реализации.

Рис. 4. Идентификация режима неустойчивости сигнала мощности РУ методом фрактального анализа в эксперименте 2

Локальные фрактальные особенности сигнала важны для временного определения (прогнозирования) начала/окончания особенностей в состоянии физической системы, слежения за ней в реальном временном режиме, использования как дополнительного средства диагностирования, независимого от человеческого фактора.

Приложение

МОДЕЛИРОВАНИЕ ФРАКТАЛЬНОГО БРОУНОВСКОГО ДВИЖЕНИЯ

Смоделируем в MS Excel одномерное фрактальное (с параметром Херста H, где 0 < H < 1) броуновское движение (ФБД) частицы с координатным положением x (t) для текущего времени t.

Используем следующее задание данных и характеристики [4] броуновского движения:

1. Начальное положение x (10) = 0;

2. Приращения Ax = x (t2) - x (t1) являются нормальными случайными величинами X с математическим ожиданием <(Ax)> = 0 и дисперсией <(Ax)2> = s2| A t|2H для любых 11, 12 с заданным приращением At = t2 - t1 и константой а (масштабным параметром значений, а2 - дисперсия приращений Ax при At = 1).

Фрактальное броуновское движение с параметром H = 0.5 совпадает с классическим броуновским движением (стохастическим винеровским процессом), для которого моделирование дискретного ряда может быть осуществлено по рекуррентной формуле:

x (t,) = x (tM) + N (0, 1), i = 1, 2, ..., n, где N(0, 1) - нормально распределенная случайная величина с нулевым средним и единичной дисперсией, в Excel моделируется функцией НОРМОБР(СЛЧИС(); 0; 1). Локальный вид модельного классического броуновского движения представлен на рис. П-1.

Расчет индекса вариации для окна размером в 1025 точек определил значение v = 0.502. Также параметр Н может быть найден из закона дисперсии, вычисляя логарифм от СКО приращений Ax для различных величин интервалов A t.

log(CKO(D x)) = const + H • log|A t|

Определяя угловой коэффициент тренда кривой по точкам {log|A t|, log(CKO(A x))} при | A11 = 1, 2, ..., 10, для модельного классического винеровского ряда получим приближение Н = 0.499.

Моделирование профиля обобщенного ФБД на основе RMD-метода (случайного перемещения средней точки) [4] для заданного параметра Херста H = 0.9 на диадических рациональных числах вида k/2n временного интервала [0, 1] представлено графиком на рис. П-2.

Расчетное значение индекса вариации для окна размером в 1025 точек v = 0.094. Закон дисперсии для данного метода моделирования выполняется лишь частично, что привносит отклонение в вычислении показателя H = 0.862 с ошибкой 4 %.

Рис. П-1. Винеровский случайный процесс (параметр H = 0.5)

0 Ф, 0.2 0,-5 0.4 0.5 0.6 0.7 0.3 0,9 1

\

\

N {

Ч /

Рис. П-2. Вид графика ФБД (параметр H = 0.9)

Литература

1. Калинушкин А.Е., Митин В.И., Семченков Ю.М. и др. Опыт создания и внедрения современной системы внутриреакторного контроля (СВРК-М) для реакторов ВВЭР-1000. Тезисы к докладу на МНТК-2010.

2. АркадовГ.В., Павелко В.И., Финкель Б.М. Системы диагностирования ВВЭР. М.: Энергоато-миздат, 2010.

3. СемченковЮ.М., МильтоВ.А., ПинегинА.А., ШумскийБ.Е. Анализ шумов нейтронного потока, вызванных флюктуациями параметров теплоносителя в активной зоне ВВЭР. Атомная энергия. Т. 103, вып. 5, 2007. С. 283-286.

4. Кроновер Р.М. Фракталы и хаос в динамических системах. Основы теории: Пер. с англ. М.: Постмаркет, 2000.

5. Дубовиков М.М., Старченко Н.В. Экофизика и фрактальный анализ финансовых временных рядов. УФН. Т. 181, № 7, 2011, С. 779-786.

6. Федер Е. Фракталы. М.: Мир, 1991.

7. Смоленцев Н.К. Основы теории вейвлетов. Вейвлеты в MATLAB. М.: ДМК Пресс, 2005.

Abstract

The paper describes the use of fractal geometry characteristics for analysis of experimental neutron data which is recorded by measuring system during testing of reactor plant core. Moments of qualitative changes in the system are determined by using the wavelet filtration of noise interferences in the measuring channel and calculation of the local fractal index over relatively short discrete sequence of equidistant time signal. The paper presents results for reactor behavior with thermal-hydraulic instability of coolant circulation and related stochastic oscillations of neutron power signals.

Keywords: Fractal dimension, index of variation, instability, time series, signal processing, wavelet transform.

Паньгин Александр Викторович, инженер Центра информационных технологий, г. Сосновыый Бор, pang@sbor.net

© Наши авторы, 2012. Our authors, 2012.