Идентификация характеристик стохастических систем на основе методов регуляризации и сингулярного анализа

Во введении обосновывается актуальность темы диссертационной работы, формулируются цель и задачи исследования, определяются новизна полученных результатов и их практическая ценность для технических и медицинских приложений, излагаются положения, выносимые на защиту.
В первой главе рассматривается современное состояние стохастического анализа и идентификации законов распределения случайных величин и процессов. Отмечается, что многим сложным структурам и процессам соответствует полимодальная характеристика распределений. Актуально дальнейшее развитие методов и алгоритмов идентификации систем с полимодальными характеристиками для создания алгоритмической и программной базы диагностических и экспертных систем, критичных к высокому разрешению характеристик и времени принятия решений.
Во второй главе анализируется и модифицируется базовый метод идентификации плотностей вероятностей случайных величин с тригонометрическим функциональным базисом. Исследуются ограничения на метод. Существующий базовый метод устойчив, если число обусловленности решаемой системы линейных уравнений не превышает 103÷104.
Для устойчивой идентификации коэффициентов разложения искомой
плотности по выбранному базису можно использовать классический метод
регуляризации А.Н. Тихонова. Однако для минимизации сглаживающего
функционала А.Н. Тихонова нужно априорно задать в некоторой норме
погрешность  правой части и погрешность матричного оператора h . При
этом выбор параметра регуляризации  осуществляют в соответствии с
принципом обобщенной невязки как решение нелинейного уравнения 22
()|| y || (h||||) 0 в итерационном цикле, согласовывая с Rh 
невязкой заданные погрешности. Здесь  , y , ,    , h – исходные
R
данные, Rh , y – приближенно заданные матрица и правая часть системы
уравнений;  – решение нормальной системы при конкретном значении  . Таким образом, для метода А.Н. Тихонова принципиально задание
совокупности    , h. Оптимальный параметр регуляризации  находят,
например, методом хорд по рекуррентной формуле с заданием начального значения α0.
Для решения СЛАУ на этапе нахождения коэффициентов разложения плотности вероятности также подходит метод псевдообращения матриц на основе сингулярного разложения.
Псевдорешение системы уравнений R  y USVT  находят по формуле opt  R y V SUT y , где R – обратная матрица Мура-Пенроуза; U,V –
h

ортогональные матрицы; S  – матрица «псевдообратная» к матрице сингулярных чисел S .
Псевдообращение – эффективный алгоритм, но в случае, когда в спектре распределения сингулярных чисел наблюдается быстрый спад к нулю (свидетельство близости к вырожденности и заведомо плохой обусловленности системы) вычислительные погрешности в значениях малых сингулярных чисел могут приводить к большим искажениям оптимального решения.
Поэтому в диссертации разработан подход, в значительной степени снимающий проблемы обоих методов, а именно: выбран алгоритм ЕС- регуляризации, не требующий априорного задания погрешностей оператора и правой части (что затруднено в методе идентификации); во-вторых, разработана эффективная регуляризация сингулярного разложения. В методе ЕС-регуляризации (эпсилон-структуризация) решается полная проблема собственных значений, а далее на  – сети вариацией  и шага сети, находится приближенное устойчивое решение. При этом автоматически определяется оптимальная величина параметра регуляризации  .
Оба алгоритма применены для идентификации законов распределения и позволили решать задачи с числом обусловленности систем в диапазоне 104÷1017.
Таким образом, нужно найти решение системы уравнений R φ ≈ y, где R
– соответствующая матрица размера l×n, y – вектор, составленный из l
независимых замеров наблюдаемой функции в условиях ограничения
2n
  2  r2 , что эквивалентно минимизации сглаживающего функционала
i1 i
M() yR22,
где  – соответствующий ограничению параметр регуляризации (множитель Лагранжа). Вектор  здесь представляет искомые коэффициенты разложения для идентифицируемой плотности распределения.
В терминах матричных операторов решается регуляризованная нормальная система
(RTRE)RT y.
Алгоритм ЕС состоит из ряда шагов.
(1)
Вначале вводится  – сеть в пространстве решений системы, состоящая из векторов следующего вида
j Vn,
j1 
j j
где  j и  j – собственные числа и собственные векторы матрицы RTR,  j –
произвольные целые числа,  – скалярный параметр, задающий шаг  – сети. При фиксированном значении параметра  решается нормальная система (1), ищется узел  – сети, ближайший к найденному решению V .
Далее оценивается качество решения  при фиксированном  , которое складывается из двух составляющих: оценка величины функционала
среднего риска, достигнутая на векторе V  плюс оценка изменения величины этого же функционала при переходе от  к V . Минимизацией функционала
по  и  определяется величина параметра регуляризации  и величина ограничения r, при которой полученное решение  оптимально при
заданном объеме экспериментальных данных.
В предложенном варианте модификации найденное оптимальное
значение  будет считаться масштабным множителем для задания равномерно распределенной случайной величины на отрезке [0,1]. В результате получается упорядоченный по убыванию компонент случайный вектор параметров регуляризации vect , максимальная проекция которого
приблизительно равна (но не превосходит)  . Далее этот вектор формирует диагональ квадратной диагональной матрицы, которая подставляется в (1) вместо E для нахождения решения  и производится вычисление величины эмпирического риска на полученном решении
IE 1 yR 2. (2) l
Рассмотренный выше алгоритм регуляризации обобщен на случай нормального распределения локальных параметров регуляризации. Алгоритм сохраняет этап формирования локальных параметров регуляризации и содержит процедуру отсечения слишком малых значений.
Обобщенный ЕС-алгоритм верифицирован на тестовых задачах с матрицами повышенной размерности. Приведем решение интегрального уравнения Фредгольма I рода по искаженным случайной помехой правой части. Интегральное уравнение сводится к алгебраической системе с числом уравнений M = 210 и числом неизвестных K = 126. Результаты решения уравнения показаны на рисунке 1.
а) б) в)
Рисунок 1. Влияние регуляризации на решение плохо обусловленной системы
а) аналитическое решение интегрального уравнения; б) SVD-решение без регуляризации; в) решение алгоритмом ЕС-регуляризации
Помеха в правой части имеет нормальный закон распределения с нулевым математическим ожиданием и с.к.о. равным 0,001. Элементы матрицы вычислялись по формуле:
11

A (I,J) = [0.164. (I-1) – 0.328. (J-1)]+
I = 1,…, 210; J = 1,…,126; […]+ = {x, x ≥ 0; 0, x < 0}. Матрица A предельно плохо обусловлена: cond_A = Inf (бесконечно большое число в MATLAB), rank(A) = 105. Время вычисления регуляризированного решения не более 1.5 с. Для идентификации законов распределения характеристик мембранных структур и реализаций ЭГЭГ-сигнала (глава 4), в основном, использовался базис тригонометрических функций. Развитием метода является включение в алгоритм идентификации функционального базиса в виде ортогональных сферических полиномов Лежандра. Исследования на тестовых примерах показывают: полиномы Лежандра позволяют более эффективно восстанавливать полимодальные законы распределения с негладкими плотностями. На рисунке 2 показаны результаты идентификации случайной величины с бимодальной равномерной плотностью распределения для тригонометрического базиса и базиса Лежандра. Случайная величина моделирует множество Кантора первого порядка. а) б) Рисунок 2. Результаты идентифиции плотности бимодального равномерного распределения с помощью: а) тригонометрического базиса; б) полиномов Лежандра. Объем выборки L = 1200 отсчетов В базисе полиномов Лежандра (график б) можно получить более крутые фронты распределения и меньшие выбросы вершинных участков, чем в тригонометрическом базисе (а). Оптимальное число членов разложения p(x) в обоих случаях примерно равно (15 и 17). Идентификация по гармоническому базису имеет небольшое преимущество в скорости вычислений. На рисунке 3 показаны результаты восстановления плотностей нормального (с малой дисперсией σ = 0.05) и распределения с тремя выраженными модами. Точность идентификации гладких распределений у обоих базисов одинакова. Рассмотренные функциональные базисы ортогональны, имеют высокое разрешение и могут эффективно идентифицировать полимодальные распределения по выборкам как малого (50-100 отсчетов), так и большого объема (до 10000 отсчетов). 12 а) б) Рисунок 3. Графики идентифицированных тестовых распределений с гладкими плотностями: a) нормальное распределение; b) полимодальное распределение. Объем выборки L = 1200 отсчетов. Синяя линия: базисные функции - нормированные полиномы Лежандра. Зеленая – тригонометрический базис Вместе с тем следует отметить, что тригонометрический базис, обладая высоким быстродействием, имеет недостаток, выражающийся в том, для распределений с равномерной плотностью или уплощенного вида дает трудно устранимую колебательную составляющую (рисунок 4). а) б) Рисунок 4. Графики идентифицированных плотностей равномерного распределения на отрезке [0,1] (объем выборок соответственно 1200 и 4800); жирная синяя линия – базисные функции – полиномы Лежандра, тонкая красная – тригонометрический базис С восстановлением таких плотностей хорошо справляется метод на основе базиса Лежандра. Например, для равномерного распределения идентифицированный энтропийный коэффициент равен 1,732 (при теоретическом значении 1,730), математическое ожидание равно 0,510. В диссертации исследован компромиссный вариант, сочетающий преимущества обоих базисов. В качестве базиса разложения плотности вероятности предложен простой степенной полином. Указанный базис не является ортогональным, восстановление коэффициентов полинома становится потенциально неустойчивой, некорректной процедурой. Для решения этой задачи использован регуляризирующий ЕС-алгоритм, рассмотренный и модифицированный в разделах 2.3 и 2.4. 13 Вариант исследования распределения с двумя разнесенными модами (число отсчетов L = 2400); показан на рисунке 5. В этом примере алгоритмически объединены нормальные распределения с центрами группирования в точках x = 0 и x = 12. Встроенный в программу идентификации регуляризирующий модуль позволяет восстановить обе моды, хотя и не в полной мере. Вычислительный алгоритм работает в «жестких» условиях больших значений чисел обусловленности и практически абсолютной вырожденности решаемых СЛАУ. а) б) Рисунок 5. Идентифицированные плотности распределения бимодальной случайной величины с модами, распределенных по нормальному закону; а) функциональный базис – полином 12 степени; б) тестовое распределение Приведем «эволюцию» значений чисел обусловленности последовательно решаемых систем: 11.721; 1451.7; 75989; 1.4709e+07; 2.3154e+08; 1.5716e+10 . . . 1.7857e+12. Очевидно, что достаточно успешно решена существенно плохо обусловленная задача. Найденный оптимальный параметр регуляризации αопт = 4.40e-04, время идентификации составляет примерно 3 с. Вывод: степенной базис имеет высокое быстродействие и эффективен при восстановлении одномодальных уплощенных кривых распределения (равномерных, типа Шапо, трапецеидальных). Регуляризация обеспечивает устойчивость при полимодальной идентификации. В третьей главе предложен способ регуляризации SVD-разложения матриц, позволяющий сделать устойчивым решения систем с большой размерностью (несколько сот уравнений и неизвестных) и сложной структурой сингулярного спектра (быстрый спад сингулярных компонент к нулю). Известный алгоритм регуляризации состоит в замене матрицы R в системе уравнений R USVT   y на приближенную Rp U SpVT , где Sp diag(1,2,...,l,l,l,...,l). Предложено неизвестный малый пороговый параметр l находить по алгоритму ЕС-регуляризации, когда значение l 14 «уподобляется» параметру регуляризации  . Исходная система R  y приводится к нормальному виду RTp Rp  RTp y (3) при этом, в отличии от классического метода А.Н. Тихонова, информация о характеристиках регуляризации скрыта внутри матрицы Rp . Эта матрица аппроксимирует сингулярное разложение матрицы R : R p  U S p V T , где Sp  diag(1,2,...,,,,...,). В предложенном способе выбор оптимального  выполняется с использованием ЕС-алгоритма. Алгоритм регуляризации состоит из следующих этапов. Сначала вычисляется «строгое» сингулярное разложение матрицы R USVT , где S содержит изначальный спектр с малыми сингулярными значениями. Затем следуют шаги: Шаг 1 – без изменений для исходной матрицы RT R (задается  - сеть в пространстве решений системы); Шаг 2 - при фиксированном значении параметра  решается нормальная система (3) и ищется узел  - сети, ближайший к найденному решению  ; Шаг 3 – оценивается качество решения  при фиксированном  . Шаг 4 - уменьшается значение  , формируется новое приближение к матрицам S и R . Определяется число малых элементов  i равных  в S p . Элементы матрицы S , меньшие или равные  , заменяются на значение  , вычисляется новая матрица Rp , затем происходит переход ко второму шагу. Шаги (2 – 4) повторяются до минимизации функционала среднего риска. Найденное регуляризированное SVD-разложение используется в алгоритме псевдообращения для нахождения устойчивых приближенных решений систем, если норма вариации правой части y не превосходит допустимого значения. Для контроля и прогнозирования устойчивости решений систем уравнений развит новый подход, заключающийся в нахождении законов распределения самих сингулярных чисел. Если в цикле случайных возмущений матрицы R идентифицируется распределение сингулярного числа  i с несколькими значимыми модами, то семейство R*  y с наиболее вероятными псевдообратными матрицами T может порождать более одного вероятного псевдорешения и множественность центров группирования. В силу убывания диагональных элементов матрицы S , псевдообратная матрица R будет более зависеть от тех элементов матрицы S, которые имеют меньшие значения. Идентифицированные плотности распределения малых сингулярных значений служат для контроля и прогнозирования устойчивости нормальных псевдорешений систем уравнений.   V  RSU В диссертации стохастический анализ сингулярных чисел матрицы основывается на методе Монте-Карло и идентификации их плотностей вероятности. Рисунок 6 иллюстрирует распределения сингулярных чисел тестовой матрицы A в условиях гауссова возмущения квазидиагонали по формуле aii + 0.5*randn. Ранг матрицы равен 2, число обусловленности: 3.4e+16. Видно, что плотности распределения малых сингулярных чисел имеют полимодальный характер. Матрица A и вектор правой части b имеют вид 1 2 1 3  2121 A = 3 3 3 4 , 1 1 1 2 4 1 4 3  T b = [1; 4; 5; 3; 10] . (4) SVD-разложение дает нормальное псевдорешение совместной системы: x = [ 1.101; -0.014; 1.101; -0.391 ]T. Вектор сингулярных чисел в отсутствии возмущения вычислен программой SVD и равен (8.307 6.928 0.000 0.000)T. в) σ3 г) σ4 Рисунок 6. Идентифицированные плотности распределения малых сингулярных чисел при решении возмущенной СЛАУ; объем выборки L = 400 Идентификация законов распределения сингулярного спектра выполнена методами, изложенными в главе 2. Математическое ожидание σ3 ≈ 0.38, а «средний» ранг матрицы увеличился до 4. Сильная трансформация имеет место для рассмотренной СЛАУ когда b = [1; 3; 5; 3; 10]T и система становится несовместной. В четвертой главе приведены результаты прикладных исследований и стохастического анализа топологий мембранных структур (атомно-силовая микроскопия), в гастроэнтерологии (электрогастроэнтерографический сигнал) и при верификации решений интегральных уравнений (ИУ) Фредгольма I рода, представляющих актуальные технические приложения. Отмечается, что ЭГЭГ-сигнал можно отнести к классу фрактальных случайных сигналов (самоаффинный фрактал), по другой терминологии – к 16 процессам с дробным броуновским движением. Идентификация выборок различной длины выявила полимодальность закона распределения как на отрезках, так и сигнала в целом (рисунок 7). Идентифицированные полимодальные характеристики содержат информацию о динамике состояний органов жедудочно-кишечного тракта (ЖКТ) пациента. Рисунок 7. Отрезок фрактальной реализации из 800 отсчетов и его плотность распределения В отличие от известных, предложенный подход восстанавливает эволюцию плотностей распределения практически в режиме реального времени и позволяет эффективно использовать методы энтропийного анализа. Высокая чувствительность энтропийных алгоритмов к изменению параметров биоэлектрических нестационарных сигналов делает востребованным предложенный подход для эффективного подбора лекарственных средств и индивидуальной терапии высокого качества. В результате выполненных исследований создан блок алгоритмического и программного обеспечения для создания перспективной экспертной системы персональной медицины. Полученные данные предполагается использовать также для экстренной оперативной диагностики. Результаты идентификаций отрезков на различных временных масштабах планируется применить при создании моделей, отображающих поведение системы ЖКТ человека на уровне уравнений диффузии или Фоккера-Планка с дробными производными. ЭГЭГ-сигнала Кроме того, разработанная модификация SSA-метода для векторов вложения траекторных матриц позволяет выделять опорные компоненты, сокращать порядок траекторной матрицы и сжимать фрактальный ЭГЭГ- сигнал без потерь значимой информации. В настоящее время актуальность в создании таких мониторинговых экспертных систем наблюдается в сфере космической медицины, в частности, в рамках эксперимента «Спланх», который проводит ГНЦ РФ– ИМБП РАН с помощью приборов «Гастроскан ГЭМ» или «Спланхограф» для исследования структурно-функционального состояния пищеварительной системы в условиях космического полета (проведение ЭГЭГ-диагностики в условия невесомости на МКС). В таблице 1 приведены результаты идентификации векторов вложения матрицы A(230,271) для временного ряда из 500 отсчетов ЭГЭГ-сигнала. Таблица 1 – Идентифицированные свойства векторов вложения траекторной матрицы Ганкеля Nopt Hp mx k tвыч,с Число мод в распределении 1 14 8.77 30 12 8.75 60 5 8.61 90 5 8.76 120 5 8.83 150 5 8.89 180 5 8.95 210 7 8.92 240 2 9.11 270 2 9.12 -3.6144e+003 -4.3390e+003 -4.7160e+003 -4.7343e+003 -4.8097e+003 -4.9793e+003 -5.0813e+003 -4.8936e+003 -5.0562e+003 -5.1342e+003 1.168 1.062 5 1.392 0.735 5 2.008 0.782 1 1.800 0.640 2 1.869 0.672 1 1.762 0.656 2 1.807 0.593 2 1.731 0.750 2 1.970 0.797 1 1.963 0.828 1 Здесь Nopt - оптимальное число тригонометрических функций, достигнутое при минимизации функционала; Hp - энтропия в форме Шеннона; mx - математическое ожидание; k - энтропийный коэффициент; tвыч, с - время вычисления плотности распределения. В приведенном примере алгоритм позволил сжать отрезок ЭГЭГ-сигнала в два раза практически без потери информации. Рисунок 8 иллюстрирует результаты идентификации плотностей распределения отмеченных векторов вложения матрицы Ганкеля. а) б) Рисунок 8. Плотности распределения векторов вложения матрицы A: а) K = 1, 30, 60, 90, 120; б) K = 150, 180, 210, 240, 270. K - номер столбца траекторной матрицы. Последовательность плотностей: синий, зеленый, красный, бирюзовый, розовый 18 Вектора с почти совпадающими кривыми распределения (контроль по Hp) можно опустить, сделав повторную «ганкелизацию» матрицы. В разделе 4.2 выполнена стохастическая оценка распределений пор и глубинных профилей полимерных мембран. Измерения проведены с использованием сканирующего зондового микроскопа SPM-9700 (Shimadzu, Япония) с максимальным полем сканирования размером 30×30 мкм2 в научном центре кафедры «Нанотехнологии и биотехнологии» НГТУ им. Р.Е. Алексеева (измерения по методу атомно-силовой микроскопии). Восстановленые плотности распределения размеров и глубин пор также носят полимодальный характер с числом мод от 1 до 6. Исследованы распределения размеров глубин пор для двух различных сечений исследуемой мембранной структуры ХТЗ-б-АН после газоразделения. Для конкретных технологических задач обнаруженные полимодальные характеристики можно соотнести с селективными, прочностными и другими физико-химическими свойствами синтезированных мембран. Предложенный подход позволяет выявлять влияние действующих технологических факторов и их вклад в топологические параметры и погрешность изготовления пленок. На рисунке 9 приведены графики глубинных профилей почти периодической мембранной структуры и ее идентифицированная плотность распределения при выборке в 512 отсчетов. а) б) Рисунок 9. График глубинных профилей мембранной структуры (а) и ее идентифицированая плотность распределения (б) Две основные моды значительно разнесены, а число базисных гармоник, потребовавшихся алгоритму идентификации равно 18. Время вычисления – 0,57 с. Расчет дифференциальной энтропии по идентифицированным распределениям позволил определить степень периодичности и качество мембраны по всем сечениям синтезированной структуры (256 линий сканирования). В третьем разделе главы 4 для исследования применимости ЕС- алгоритма регуляризации для важных прикладных задач решены пять уравнений Фредгольма I рода. Такие задачи возникают при математической обработке данных физических экспериментов при сейсморазведке полезных 19 ископаемых, в томографии высокого разрешения, при оптимальной закалке образцов, синтезе многослойных оптических структур. Рассмотрим одно из решенных автором интегральных уравнений (ИУ) Фредгольма I рода с гладким симметричным ядром: 11 Az(ts)21z(s)dsu(t), t[1,1] 1 Пусть z(s) задана выражением: z(s)  3.25, s  S1  [ 0.54, 0.54]; z(s)0.25, sS1 После вычисления u(t) на заданной сетке значений, полагаем z(s) неизвестным и находим решение ИУ методом регуляризации. При этом количество узлов сетки по s и t равно соответственно 71 и 91. Число обусловленности матрицы аппроксимирующей системы алгебраических уравнений cond(A) = 3.6e+017, ранг = 25. Достигнутое ЕС- алгоритмом оптимальное значение параметра регуляризации в скалярном варианте равно 1.86e-008. Вычисления проведены в двух вариантах: а) присутствуют только ошибки аппроксимации и округления; б) моделируются ошибки в исходных данных для правых частей ИУ. Возмущение правой части выполнялось для каждой проекции по нормальному закону с нулевым математическим ожиданием и с.к.о. sigma = 0.0003. Кроме плохой обусловленности (числа обусловленности более 1017), в приведенных примерах имеет место быстрое убывание значений сингулярного спектра. Результаты решения приведены на рисунке 10. а) б) в) г) д) е) Рисунок 10. Исходные данные и результаты решения ИУ: (а) точное решение; (б) точная правая часть; (в) решение ИУ SVD-алгоритмом; (г) решение скалярным ЕС-алгоритмом; (д) решение векторным ЕС-алгоритмом; (е) решение SVD-алгоритмом при возмущении правой части 20 Для сравнения показаны графики решений по методу SVD-разложения без регуляризации. SVD-алгоритм обеспечивает более крутые фронты, но приводит к пилообразным искажениям z(s) при малейших погрешностях в правой части (график-е). ЕС-алгоритм дает устойчивые, но более сглаженные решения. Время вычисления – не более 2 с. В работе приведено также решение интегрального уравнения в задаче измерения параметров ускоряемого пучка заряженных частиц в ускорителе с помощью второго пучка частиц, распределение плотности частиц в котором известно (данные Института ядерной физики им. Г.И. Будкера СО РАН). Алгоритм ЕС- регуляризации позволяет найти устойчивое решение ИУ и восстановить неизвестную стохастическую характеристику (плотность распределения ускоряемого пучка) в сложном и дорогостоящем эксперименте. В заключении приведены основные результаты работы: 1. Показана актуальность исследований для стохастических явлений и структур с полимодальными, негауссовыми законами распределения. 2. Разработан алгоритм регуляризации для идентификации плотностей распределения полимодальных характеристик на основе модификации алгоритма ЕС-регуляризации. 3. Предложен алгоритм идентификации полимодальных распределений с функциональным базисом в виде ортогональных сферических полиномов Лежандра. Выполнен сравнительный анализ устойчивости и точности идентификации трех базисов: тригонометрического, на основе многочленов Лежандра и простого степенного базиса. 4. Модифицированный алгоритм протестирован при решении плохо обусловленных систем повышенной размерности. 5. Предложена регуляризация SVD-разложения с глобальным параметром регуляризации и реализована его верификация. 6. Выполнена идентификация плотностей распределения сингулярных значений в методе возмущения для плохо обусловленных матриц общего вида. 7. Обобщен SSA-метод (Singular Spectrum Analysis) для исследования траекторных матриц случайных процессов. 8. Метод ЕС-регуляризации и его модификация доказали высокую эффективность при решении сложных интегральных уравнений Фредгольма первого рода в прикладных технических задачах. 9. Результаты исследований в технических и медицинских приложениях подтверждают работоспособность предложенных алгоритмов и методов.