Разработка алгоритма идентификации

Разработка алгоритма идентификации

СОДЕРЖАНИЕ

Определения, обозначения и сокращения

Введение

1. Научная аппаратура СОКП и исходные данные для анализа

1.1 Состав и принцип работы аппаратуры

1.2 Исходные научные данные для анализа

2. Выбор параметров идентификации

2.1 Выбор параметров корреляционного анализа

2.2 Выбор параметров Фурье-анализа

3. Алгоритм корреляционного анализа

3.1 Разработка алгоритма корреляционного анализа

3.2 Применение алгоритма корреляционного анализа

Заключение

Список использованных источников

Приложение 1. Реализация алгоритма Фурье-анализа на языке С++

Приложение 2. Реализация алгоритма корреляционного анализа на языке С#

Определения, обозначения и сокращения

АСУ - ассенизационное санитарное устройство

АР - автономный регистратор

БВК - блок вакуумных клапанов

БПАС - блок преобразования акустических сигналов

ВГК - внешний гидравлический контур

ГЖТ - газожидкостный теплообменник

ИУС - информационно-управляющая система

КА - космический аппарат

КС - комплексный стенд

КЭ - космический эксперимент

МКС - Международная Космическая Станция

НА - научная аппаратура

НПО ИТ - Научно-производственное объединение измерительной техники

ПО - программное обеспечение

ПрК - переходная камера

ПхО - переходный отсек

РКК - Ракетно-космическая корпорация

РО - рабочий объём

РС - российский сегмент

СКВ - система кондиционирования воздуха

СМ - служебный модуль

СОКП - средства оперативного определения координат точки пробоя

СОТР - система обеспечения теплового режима

СПН - сменная панель насосов

СУБА - система управления бортовой аппаратурой

ЦПК - Центр подготовки космонавтов им. Ю.А. Гагарина

ЦУП - центр управления полётом

УТМ - учебно-тренировочный макет

ФГУП ЦНИИмаш - Федеральное государственное унитарное предприятие Центральный научно-исследовательский институт машиностроения

Введение

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

Контроль за космическим мусором вокруг МКС идет постоянно, и при необходимости специалисты ЦУП «заставляют» МКС делать маневр «Уклонение от осколков».

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

Печальный опыт разгерметизации модулей станции в истории пилотируемой космонавтики был, когда произошла разгерметизация модуля «Спектр» космической станции «Мир». Тогда найти место разгерметизации модуля и восстановить его работу не удалось.

В связи с этим при проектировании новых пилотируемых космических модулей и эксплуатации уже существующих должна учитываться возможность пробоя гермооболочки в результате столкновения с техногенной частицей или микрометеороидом. Вероятность непробоя гермооболочки МКС даже при усиленной экранной защите модулей оценивается величиной на уровне 0,85 в течение 15 лет эксплуатации.

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

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

В настоящее время на РС МКС проводятся космические эксперименты «Отклик» и «Пробой». В КЭ «Отклик» исследуется возможность определения места соударения гермооболочки космического модуля путём анализа волн вибраций, распространяющихся по самому гермокорпусу модуля, а в КЭ «Пробой» - посредством анализа акустических волн, распространяющихся в воздушной среде модуля.

Космический эксперимент «Пробой» проводится с целью верификации метода оперативного определения координат точки пробоя. Подзадачами КЭ «Пробой» являются изучения собственной акустической обстановки и идентификация постоянно или периодически действующих на борту РС МКС источников импульсного акустического шума. Данная идентификация необходима в связи с потребностью отсеивать известные сигналы работы бортовой аппаратуры от потенциально возможных сигналов реального пробоя. Изучение акустической среды в РС МКС необходимо для учёта всех параметров при проектировании штатной системы, максимизации точности и устойчивости решений.

Шумящее оборудование на борту РС МКС включает в себя насосы СПН и ГЖТ ВГК СОТР, клапаны БВК, элементы системы АСУ, циркуляционные и комфортные вентиляторы, компрессор и вентилятор СКВ, различная научная аппаратура и многое другое. Также дополнительными источниками акустических импульсов являются действия самих космонавтов: занятия на беговой дорожке, приготовление пищи, работа с приборами, открытие и интерьерных панелей и дверей кают, разговоры с Землёй и между собой. Спектры звуков защёлок интерьерных панелей были отдельно исследованы на УТМ в ЦПК им. А.Ю. Гагарина.

Распознавание сигналов широко используется в различных сферах человеческой деятельности: в сейсмологии [3, 4], океанологии [5], кардиологии, сфере IT [6, 7, 8] и многих других.

При цифровом распознавании сигналов используются различные методы: сингулярные разложения [9], метод главных компонент [10], спектральный анализ и анализ фурье-методами [11, 12], корреляционные методы [12], вейвлет-преобразования [13, 14].

Целью данной работы является идентификация источников импульсного акустического шума, действующих на борту РС МКС, для дальнейшей доработки научной аппаратуры «СОКП» и создания штатной системы.

Задачей работы является разработка алгоритма, позволяющего эффективно различать импульсные шумы от различной бортовой аппаратуры, его реализация в программном коде и применение для анализа научной информации, полученной в ходе КЭ «Пробой».

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

В первом разделе работы описан научная аппаратура «СОКП», её устройство и принцип работы, а также формируемые ею в ходе КЭ «Пробой» данные для последующего анализа.

Во втором разделе проанализированы критерии, по которым следует идентифицировать источники импульсного акустического шума, произведён выбор параметров анализа.

В третьем разделе описан алгоритм, который позволяет разделять импульсные сигналы на группы, различного происхождения, и приведены результаты его применения к научной информации, полученной в ходе КЭ «Пробой» в период с декабря 2014 года по ноябрь 2015 года.

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

1. Научная аппаратура СОКП и исходные данные для анализа

.1 Состав и принцип работы аппаратуры

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

Результаты космического эксперимента планируется использовать при создании штатной системы оперативного определения координат места пробоя.

Принцип действия системы оперативного определения координат пробоя основан на измерении разновременности прихода акустической волны от точки удара к акустическим датчикам (микрофонам), расположенным в различных точках гермоотсека по определенной схеме (на способ и схему размещения датчиков имеются патенты [1, 2]), и расчете координат источника акустической волны на основании этих данных.

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

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

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

Метод определения координат пробоя был успешно подтвержден в экспериментах, проведенных на КС СМ в РКК «Энергия» и на УТМ СМ в ЦПК им. Ю.А. Гагарина.

Постановщиком эксперимента «Пробой» и разработчиком научной аппаратуры является ФГУП ЦНИИмаш. В разработке и изготовлении имитатора пробоя принимало участие ОАО «РКК «Энергия», а в изготовлении вторичных источников питания - ОАО «НПО ИТ».

Научная аппаратура КЭ «Пробой» состоит из двух независимых регистрирующих блоков: БПАС (рисунок 1), АР (рисунок 2) с вторичными источниками питания, комплекта микрофонов в адаптерах, имитатора пробоя (рисунок 3).

Рис. 1. БПАС

Рис 2. АР

Рис. 3. Имитатор пробоя

БПАС размещается в РО СМ, АР используется как переносной регистратор и размещается поочередно в ПхО СМ или ПрК СМ.

К каждому из регистрирующих блоков подключено по 6 датчиков. В качестве чувствительных элементов используются фазированные ICP-микрофоны массива типа 4958 фирмы “Bruel&Kjaer”. Микрофоны закреплены в специальных адаптерах, которые защищают их от случайных повреждений. Адаптеры крепятся на панелях интерьера или вблизи гермооболочки с помощью велкро (рисунок 4).

Рис. 4. Микрофон в сборе закреплен на панели интерьера

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

На рисунке 5 представлена структурная схема НА «СОКП» взаимодействия с системами СМ.

Рис. 5. Структурная схема НА «СОКП»

Электрическая схема подключения НА «СОКП» приведена на рисунке 6.

Рис. 6. Схема подключения НА «СОКП»

Основные технические характеристики НА «СОКП» приведены в таблице 1.

Таблица 1

Технические характеристики НА «СОКП»

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

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

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

Основные затруднения при определении момента прихода переднего фронта волны от гермооболочки связаны с прохождением звука в рабочее пространство модуля через приборный отсек, отделенный от рабочей зоны обитания космонавтов панелями интерьера толщиной от 8 до 20 мм, изготовленными из вспененного полистирола, обшитого тонким слоем алюминиевой фольги толщиной ~ 0.3 мм и оклеенными ворсовой тканью типа «велкро». Наземные эксперименты показали, что при прохождении акустической волны через панель интерьера ее амплитуда снижается на 30-40 дБ.

Учитывая, что уровень звука, генерируемого при пробое, составляет не менее 140 дБ на расстоянии 1 м от точки пробоя, а уровень фоновых шумов в рабочем отсеке составляет 60-70 дБА, можно достаточно надежно фиксировать акустическую волну пробоя после ее прохождения сквозь панели интерьера.

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

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

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

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

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

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

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

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

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

Все результаты коротких и длительных сеансов верификации записываются в файлы на флэш-накопители, которые возвращаются на Землю с экипажем.

1.2 Исходные научные данные для анализа

Акустический импульс (от лат. impulsus - удар, толчок) - это бегущая звуковая волна, имеющая характер резкого кратковременного изменения давления. Такие звуковые волны создаются, например, взрывом, искровым разрядом, соударением тел. Каждый такой импульс содержит как область повышенного, так и область пониженного давления. Спектр такого акустического импульса сплошной, с максимумом в области частот, период которых близок к длительности импульса.

В 2008 году проводились исследования пробоя фрагмента гермооболочки модуля на установке ЦНИИмаш с помощью газокумулятивной пушки, в стволе которой снаряд из стали или алюминия массой до 0,2 г разгонялся до 5 км/c. На рисунках 7 и 8 показан импульс в разных масштабах времени, а на рисунке 9 - спектр данного импульса.

Рис. 7. Импульс пробоя гермооболочки высокоскоростным снарядом

Рис. 8. Развёртка импульса пробоя гермооболочки

Рис. 9. 1/3-октавный спектр импульс пробоя гермооболочки

алгоритм корреляционный фурье язык

Из рисунка 8 видно, что длительность основного импульса составляет порядка t_{импульса} = 0,01 с.

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

Первый этап КЭ «Пробой» проводился на борту СМ РС МКС в трёх экипажами трёх миссий:

миссия 41-42 с декабря 2014 по февраль 2015;

миссия 43-44 с мая 2015 по июнь 2015;

миссия 45-46 с октября 2015 по ноябрь 2015.

В каждой миссии проводилось несколько коротких сеансов верификации и длительные сеансы непрерывного включения аппаратуры СОКП. Всего суммарная длительность длительных сеансов составила 91,66 суток.

В ходе первого этапа КЭ «Пробой» в РО СМ РС МКС непрерывно находилось 6 микрофонов, подключенных к БПАС. Ещё 6 микрофонов, подключенных к АР попеременно располагались в ПхО и ПрК. Микрофоны располагались согласно запатентованной схеме [1, 2], а также указаниям специалистов по бортовой аппаратуре для обеспечения комфортного взаимодействия с приборами и системами служебного модуля. В миссии 45-46 эксперимент проводился с проводной синхронизацией, для чего на 7 и 8 каналы БПАС заводились данные с АР. Координаты микрофонов приведены в таблице 2 и на рисунках 10 и 11. На рисунках 10 и 11 также представлен внешний вид программы [3].

Таблица 2

Координаты микрофонов на борту СМ

Рис. 10. Координаты датчиков (вид СМ сверху) и внешний вид программы.

Рис. 11. Координаты микрофонов (вид на правый борт СМ сверху).

Данные с микрофонов снимались с частотой опроса f_опроса = 50 кГц синхронно по всем каналам. Из этих данных формировались буферы по N_буфера = 1024 показания по каждому каналу. Значит, каждый буфер соответствует сигналу длительностью t_буфера = 0,02048 с. Как видно, t_буфера примерно вдвое превышает длительность t_{импульса}.

В каждом буфере вычислялся суммарный уровень звукового давления. Эти суммарные уровни записывались в отдельные файлы. Циклограмма типовых рабочих суток экипажа МКС (12 февраля 2015 года) по датчику М4 представлена на рисунке 12.

Рис. 12. Зависимость уровня звука (в дБА) от времени суток (в ч)

Как видно из рисунка, максимального уровня в 110 дБА шум достигает в часы работы экипажа. В ночное же время над уровнем фона превышают в основном только импульсы от срабатывания клапанов БВК каждые 15 минут.

Если суммарный уровень превышал уровень фоновых шумов минимум на 6 дБ хотя бы по двум датчикам, этот буфер, а также один предыдущий и два последующих записывались в файл импульса. Кроме самого сигнала из 4×N_буфера = 4096 точек в файлы записывались дата время регистрации этого импульса с точностью до 1 мс. На рисунке 13 приведен вид программы LoadSignal [12] с загруженной записью импульса.

Рис. 13. Порция сигнала в программе LoadSignal

С 15:13 13 ноября 2015 года по 03:56 14 ноября 2015 года производилась регистрация импульсов по превышению на 6 дБ не по двум, а по одному каналу. Это привело к аномальному количеству записей импульсов: за t₂ = 91,66 суток длительных сеансов было зарегистрировано m₂ = 387 000 импульсов, в то время как за t₁ = 0,53 суток 13-14 ноября 2015 зарегистрировано m₁ = 725 700 импульсов. Индексы в данном случае обозначают количество датчиков, по которым происходила регистрация импульса.

Аномальность числа импульсов m₁ выражается в том, что за время t₁ суммарное время всех записанных импульсы составляет

Подставив значения, получим T₁ = 0,69 суток, что превышает t₁. Это можно объяснить тем, что фактом, что дискретность регистрации импульсов равна длительности одного t_буфера, а размер записи превосходит в 4 раза это время. Таким образом, может происходить перекрытие соседних записей.

Для углубления знаний об акустической обстановке на борту РС МКС эти избыточные данные после удаления зон перекрытий были соединены в большие wav-файлы. Таким образом, удалось собственным слухом почувствовать все события, произошедшие на борту за тот период. В частности, ночью удалось зафиксировать периодические изменения в акустической обстановке, вызываемые попеременной работой одного или двух роторов, обеспечивающих циркуляцию воздуха. Циклограммы ночных записей в 1/3-октавных спектральных полосах представлены на рисунке 14.

Рис. 14. Циклограммы ночных часов в 1/3-октавных полосах частот

Из рисунка видно, что на некоторых частотах между П-образными подъёмами, соответствующими срабатываниям клапанов БВК, уровни звукового давления в некоторых 1/3-октавных полосах частот заметно меняются. Особенно это касается 1/3-октавной полосы частот с центральной частотой 3150 Гц, в которой изменение уровня звукового давления достигает 7 дБ.

Для более подробного анализа были построены узкополосные спектры с разрешением 1 Гц двух последовательных периодов работы роторов, и произведено вычитание одного из другого. Таким образом, удалось исключить общие фоновые шумы из спектров каждого интервала и более наглядно продемонстрировать, что же именно изменилось. Результат приведён на рисунке 15.

Рис. 15. Разность спектров последовательных интервалов работы роторов БВК

Верхняя (выше 0) часть графика соответствует работе одиночного ротора, а нижняя (ниже 0) - работе пары роторов. Легко заметить, что узкие (тональные) пики, присутствующие в обоих случаях, меняют своё положение и амплитуду, в зависимости от того, сколько роторов работает одновременно. Частоты этих пиков являются кратными к собственной частоте вращения ротора. По графику можно определить, что в случае работы одного ротора частота вращения составляет 110,0 Гц, а при работе двух - 107,7 Гц. Эту разницу можно объяснить тем, что двум роторам, работающим совместно, не нужно так же много усилий, для обеспечения необходимых расходо-напорных характеристик, как одному.

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

Таким образом мы показали, что НА СОКП помимо своей основной задачи оперативного определения координат точки пробоя может выполнять функции мониторинга, показывая изменение в режимах работы шумящей аппаратуры, связанные, возможно, в том числе и с её неисправностью, а также поиске резонансов, соответствующих собственным частотам акустических объёмов, чтоб использовать эти данные при проектировании шумящей аппаратуры и элементов звукопоглащения.

Скорость звука в течение эксперимента могла меняться в зависимости от параметров окружающей среды (температуры и влажности, которые определялись согласно показаниям термогидрометра «ИВА-6» непосредственно экипажем станции и заносились в электронные протоколы эксперимента) в пределах c = 340±15 м/с. Учитывая максимальное расстояние между микрофонами БПАС М1 и М5 ∆L = 5,258 м, максимальная разница во времени прихода фронта акустического импульса определяется по формуле

и составляет 0,01546 с.

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

За время проведения длительных сеансов средняя плотность регистрации импульсов ρ

и составляет около 4 220 импульсов в сутки или 0,05 импульсов / с. Положим, что источники каждого импульса независимы (хотя это не совсем так, некоторые приборы имеют характерную циклограмму). В таком случае вероятность наложения двух импульсов определяется по формуле:

и равна 2,5∙10^-7. Учитывая, что всего было зарегистрировано m₂ = 387 000 импульсов, из них двойными были в среднем

что составляет 0,1. Таким образом, приходим к заключению, что количество наложенных импульсов практически равно нулю.

Тем самым мы приходим к следующим выводам:

)        Все импульсы, вызвавшие превышения по двум датчикам - одиночные, а значит, нет необходимости искать методы их разделения.

)        Записей сигналов импульсов в 4096 точек более чем достаточно для того, чтоб зафиксировать импульс по всем датчикам.

2. Выбор параметров идентификации

.1 Выбор параметров корреляционного анализа

Идентифицировать источник импульсного сигнала можно несколькими способами: по суммарному уровню импульса, по спектру импульса и по координатам источника импульса. Для лучшего результата будем использовать комбинацию всех этих методов.

Для сопоставления двух и более сигналов по спектру и уровню проще и нагляднее сравнивать не сами сигналы, а их спектры. Выделим в каждом сигнале спектральные полосы частот, например 1/3 октавные. Пусть есть спектры сигналов  и , где i = 1, 2, …, I (I - количество спектральных полос частот). Классически коэффициент линейной корреляции определяется по формуле:

Такой коэффициент лежит в диапазоне [-1; 1]. Случай, когда с = 1, соответствует абсолютно коррелирующим сигналам. При с = -1 сигналы анти коррелирующие. Эта формула совпадает с выражением для косинуса угла между двумя векторами  и  в многомерном ортонормированном пространстве.

Именно такая форма коэффициента линейной корреляции использовалась в работе [15]. Однако она имеет существенные недостатки. Во-первых, она учитывает только разницу в спектрах, а не в суммарных уровнях, так как на значение косинуса угла между векторами не оказывает влияние их длина. То есть нельзя различить источники с похожим спектром и разными суммарными уровнями сигнала. Во вторых, вычитание математических ожиданий  и  может сделать спектры похожими, даже если они таковыми не являлись.

Рассмотрим в таком пространстве радиус-вектора , где k = 1, 2, …, K (K - количество векторов, k - порядковый номер вектора). Концы векторов образуют некоторое множество (облако) точек в этом пространстве. Причём вектора, соответствующие спектрам сигналов от одного источника, будут указывать на близкорасположенные точки. Действительно, если источники отличаются только по уровню сигнала, но схожи по спектру, тогда радиус-вектора их импульсов будут направленны в одну сторону, но различаться по длине, и чем больше будет различие по уровню, тем дальше будут разнесены концы радиус-векторов. Если же источники близки по суммарному уровню, но различаются по спектру, тогда радиус-векторы их будут направленны в разные стороны, и чем больше будет различие по спектрам, тем дальше будут разнесены концы радиус-векторов.

Теперь передвинем начало координат в центр облака с координатами

Далее будем работать с полученным множеством радиус-векторов . Такое преобразование даст бóльшую разрешающую способность при определении косинуса угла, так как на всё том же множестве (облаке) вершин радиус-векторов будут заведомо присутствовать те, которые отличаются на углы, близкие к π. То есть, если вы стоите в центре рощи, то в какую сторону бы вы не смотрели, за вашей спиной всегда будут деревья.

Запишем выражение для коэффициентов линейной корреляции - косинусов угла между этими новыми векторами:

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

Теперь выберем множество, на котором будем сравнивать корреляционные матрицы. Возьмём в качестве первого и второго сигналов спектры от имитатора пробоя снятые в реальных условиях космического эксперимента по двум разным каналам. В качестве векторов возьмём пять спектров различных сигналов в 1/3-октавных полосах частот с центральными частотами от 400 Гц до 16 кГц.

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

В качестве четвёртого спектра возьмём так называемый «розовый шум», спектральная плотность которого обратно пропорциональна частоте. Этот выбор обусловлен тем, что в шуме бортовой аппаратуры СМ РС МКС присутствуют шумы, близкие к «розовому». На самом деле, чистого «розового шума» не наблюдается, так как «розовый шум» - это лишь математическая модель, к тому же обращающая знаменатель коэффициента линейной корреляции для . Чтоб избежать этого, мы введём в него малое возмущение - добавим 0,01 дБ в 1/3-октавной полосе частот с центральной частотой 400 Гц.

Ну и, наконец, в качестве пятого сигнала возьмём копию первого, сдвинутую на 20 дБ вниз.

Спектры всех пяти сигналов представлены на рисунке 16.

Рис. 16. Спектры тестовых сигналов

Таблица 3

Сравнение корреляционных матриц

В таблице 3 приведены спектры сигналов в числовом виде, а также корреляционные матрицы для каждого из 6 случаев. Корреляционные матрицы для удобства раскрашены в цветовые градации: чем больше коэффициент линейной корреляции (ближе к 1), тем фон зеленее, чем он меньше (ближе к -1), том фон краснее.

Из таблицы 3 видно, что методы вычисления коэффициента линейной корреляции дают высокие значения (близкие к 1) для корреляции сигналов 1 и 2, что и следовало ожидать.

Корреляции для логарифмического масштаба на множестве  не удалось разрешить ни одну группу сигналов, что, очевидно, не удовлетворяет нас. Корреляция для логарифмического масштаба на множестве  причислила сигнал 3 к группе с сигналами 1, 2, а не к группе с сигналом 4, что не удовлетворяет условиям различия источников по спектру.

Все корреляционные матрицы, построенные по линейно масштабированным множествам, причислили сигнал 5 к группе сигналов 1 и 2, не удовлетворяет условию различия источников по суммарному уровню сигнала.

Таким образом, наилучший результат показал метод вычисления коэффициента линейной корреляции по логарифмическому масштабу множеств , что и было предсказано выше. Ко всему прочему, этот метод даёт более высокое разрешение на данной выборке сигналов: значения коэффициентов линейной корреляции лежат в диапазоне от -0,8 до 1,0. В дальнейшем, будем использовать именно этот коэффициент линейной корреляции как основной.

2.2 Выбор параметров Фурье-анализа

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

Звуковая волна от источника звука распространяется в среде со скоростью звука, которая при нормальных условия в воздухе составляет около 340 м/с. Уровень звукового давления в безэховом пространстве убывает с квадратом расстояния от источника. Следовательно, передний фронт волны звукового импульса на датчики, разнесённые в пространстве, придёт в разные моменты времени и с разным суммарным уровнем звукового давления. Значит, учтя разницу в моментах времени прихода переднего фронта звукового импульса от источника и разницу в уровнях звукового давления для группы импульсов по множеству датчиков, можно определить, пришли ли импульсы от источников, локализованных в различных точках пространства или нет.

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

Вычислить время прихода фронта звуковой волны можно по записи сигнала из 4096 точек (отсчётов) так же, как это делается в программе [16] по методу, описанному в [3]. Этот метод позволяет определить, на какую именно точку (соответствующую моменту времени) пришелся момент прихода переднего фронта волны импульса. Это позволит нам корректно сравнивать уровни звукового давления на разных датчиках и одни и те же моменты времени относительно этого времени прихода распространяющегося во все стороны переднего фронта волны.

На рисунках 17 и 18 представлены графики узкополосных и 1/3-октавные спектры имитатора пробоя на фоне прочих спектров.

Рис. 17. Узкополосный спектр акустического импульса от имитатора пробоя (красным) на фоне спектров прочих импульсов.

Рис. 18. 1/3-октавный спектр акустического импульса от имитатора пробоя (красным) на фоне спектров прочих импульсов.

Для сравнения спектров и уровней сигнала были выбраны 1/3-октавные полосы частот с центральными частотами от 400 Гц до 16 кГц. Из рисунков 15 и 16, что на частотах ниже 400 Гц спектры имитатора очень близки к нижней границе чёрных графиков. То есть они тонут в фоновых шумах. В 1/3-октавной полосе с центральной частой 20 кГц наблюдается снижение уровня звукового давления по всем зарегистрированным импульсным источникам, так что эта полоса также была исключена из дальнейшего рассмотрения.

Максимальная частота (частота Найквиста) согласно теореме Котельникова для частоты опроса f_опроса = 50 кГц, не может превышать значения

что составляет 25 кГц. Минимальная же частота определяется по формуле

где N_порции - размер порции (число отсчётов), которая использовалась для быстрого преобразования Фурье. Найдём достаточный для нас размер порции N. Для большего удобства дальнейшего применения быстрого преобразования Фурье будем выбирать N среди степеней двойки.

С одной стороны, чем больше N, тем больше можно захватить нижних частот для анализа. Все прочие частоты будут кратными этой f_min, а внутри каждой 1/3-октавной полосы должна присутствовать хотя бы одна частота для суммирования. С другой стороны, N желательно делать таким, что не выходить за границу импульса, а длительность некоторых импульсы может не превышать 0,01 с.

В таблице 4 приведено сравнение порций размеров 2048, 1024, 512 и 256. Для каждой порции представлены её длительность, минимальная частота, а также количество частот, попадающих в 1/3-октавные полосы с центральными частотами от 200 до 800 Гц.

Таблица 4

Сравнение порций различного размера

Из таблицы видно, что порции в 256 точек для анализа недостаточно, так как в 1/3-октавной полосе частот с центральной частотой 500 Гц нет ни одной частоты для суммирования. Так же видно, что порция в 2048 точек является излишней, так как длительность этой порции может значительно превосходить длительность некоторых импульсов. Порция в 512 точек будем считать достаточной, так как её длительность сравнима с длительностью импульса, а все 1/3-октавные полосы частот с центральными частотами от 315 Гц и выше представлены ненулевым количеством слагаемых.

Итак, был выбран метод определения номеров отсчётов по каждому из датчиков n_j (где j - номер датчика). Было выбрано количество точек, которые нужно взять начиная с выбранного номера отсчёта, чтоб реализовать разложение импульса по быстрому преобразованию Фурье. В случае если

то есть, если n_j > 3584 записи сигнала не хватит для составления порции. В таком случае, принимаем n_j = 3584. На самом деле, такая ситуация может возникнуть только на тех датчиках, звук от акустического импульса не превысит фоновых шумов или при сбое. В штатном же режиме, как было описано ранее, при записи сигнала буфер с импульсным превышением идёт вторым из четырёх, поэтому n_j ≤ 2048.

Для каждого зарегистрированного импульса по каждому каналу (датчику) выделим таким методом порцию и преобразуем её по алгоритму быстрого преобразования Фурье, чтобы просуммировав результирующие дискретные частоты получить спектры в 1/3-октавных полосах частот.

Для выполнения этих функций в программной среде разработки LabWindows CVI 2013 была написана программа, поочерёдно применяемая к исходным файлам, содержащим записи импульсных шумов, зарегистрированных блоком БПАС в ходе прошедших длительных сеансов КЭ «Пробой». Листинг этой программы приведён в приложении 1. Среда LabWindows CVI 2013 была выбрана по двум причинам. Во-первых, расчётная часть программы [16] была реализована в этой же среде, поэтому не было проблем совместимости форматов и типов данных. Во-вторых, среда LabWindows CVI 2013 содержит встроенный математический аппарат, позволяющий автоматизировано применять такие алгоритмы, как быстрое преобразование Фурье, без подключения каких-либо внешних специализированных пакетов.

На выходе программа выдает единый файл-таблицу, в строках которой располагались даты и времена регистрации импульсов и их 1/3-октавные спектры по всем каналам. То есть получается матрица A из K строк и I×J столбцов, где K - число сигналов, I - число полос частот, J - число каналов (датчиков).

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

3. Алгоритм корреляционного анализа

.1 Разработка алгоритма корреляционного анализа

Вычислив коэффициенты линейной корреляции для каждого из K сигналов со всеми остальными по формуле:

получим квадратную симметричную матрицу C размера K×K. Симметричность её очевидна, ввиду того, что

Каждый элемент этой матрицы представляет собой число из интервала [-1; 1], причём на главной диагонали - строго единицы. В этом легко убедиться, подставив k₁ = k₂.

Для наглядности распечатаем эту матрицу в точечный рисунок таким образом, что меньшим значениям (близким к -1) элементов матрицы соответствовали более светлые точки (близкие к белым), а большим (близким к 1) - более тёмные (близкие к чёрным). Получившаяся матрица представлена на рисунке 19.

Рис. 19. Упорядоченная по датам импульсов корреляционная матрица.

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

На рисунке 20 представлено распределение значений коэффициентов линейной корреляции и этой «серой» матрице.

Рис. 20. Распределение значений коэффициентов линейной корреляции.

Из графика видно, что распределение имеет форму асимметричного колокола. Максимум приходится на значение -0,13. Значения ниже -0,92 практически не встречаются. На значении 1,00 заметен «аномальный» подъём, обусловленный наличием единичных диагональных элементов.

В идеальном случае (идеальном для корреляционного анализа, а не для жизнедеятельности экипажа), каждый прибор имел бы строго одинаковые независящие от даты и времени спектры импульсных сигналов, корреляция которых была бы равна 1, корреляция со спектрами всех остальных приборов равнялась бы 0, а все прочие значения отсутствовали. Такое распределение условно изображено на рисунке 21. В таком случае, было бы очень просто определить импульсный сигнал, как исходящий от конкретного прибора.

Рис. 21. Идеальное распределение значений коэффициентов линейной корреляции.

Это возможно, например, в том случае, когда спектры источников представляют собой узкие и высокие пики, причём частоты пиков от двух любых различных источников отличаются не менее, чем на ширину 1/3-октавной полосы частот.

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

Выберем некоторое пороговое значение коэффициента линейной корреляции c_п, такое что

Теперь мы можем путём одновременной перенумерации строк и столбцов матрицы C привести её к виду

где в белых областях значение каждого элемента не будет превосходить c_п, а вдоль главной диагонали располагаться один или несколько квадратных симметричных подматриц, таких что каждая либо состоят из единственного (и единичного) элемента либо для любых двух столбцов j и k (j ≠ k) существует такая строка i, что

В предельных случаях при c_п = 1 серая область будет представлять собой располагающиеся вдоль главной диагонали K подматриц , а при c_п = -1 - единственную подматрицу, совпадающую со всей матрицей.

Сначала возьмём c_п = -1. Затем будем постепенно увеличивать значение c_п и продолжать перенумерацию строк и столбцов, до тех пор, пока единая «серая» подматрица не «развалится» на несколько. Каждой такой отпочкованной дочерней подматрице будем присваивать значение признака c_п, при котором произошло её отделение от родительской подматрицы.

Продолжая постепенное увеличение c_п, мы будем наблюдать фрактальную картину: каждая из дочерних подматриц будет вести себя так же, как полная матрица, разбиваясь при определённых значениях c_п на ещё более мелкие подматрицы. В конце концов (при c_п = 1), когда вся матрица разобьётся на K подматриц , каждой из которых будет присвоено собственное числовое значение признака c_п.

Для реализации этого рекурсивного алгоритма в программной среде Microsoft Visual Studio 2010 была написана программа для корреляционного анализа. Листинг этой программы приведён в приложении 2. На вход программы подаётся файл с неупорядоченной (упорядоченной по времени) матрицей A, на выходе - файлы с упорядоченной матрицей C’ (в том числе и в графическом виде для наглядности) и с матрицей A’, строки которой упорядочены так же, как перенумировываются строки и столбцы матрицы C при выполнении алгоритма для получения матрицы C’. Эта матрица A’ расширена столбцом, состоящим из значений признаков, при которых происходило выделение подматриц.

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

Таким источником, например, является имитатор пробоя: от одного короткого сеанса верификации к другому менял своё местоположение согласно методике проведения космического эксперимента. Основные же бортовые источники тем и удобны в распознавании, что их местоположения постоянно относительно компоновки модуля СМ.

Другим примером могут являться парные или множественные источники импульсного акустического шума, представленные различными экземплярами одного изделия, разнесёнными в пространстве. Практика показывает, что в зависимости от конкретного изделия различные экземпляры могут иметь или не иметь значительного разброса по уровню звукового давления. В таком случае НА «СОКП» сможет улавливать эти различия или их отсутствие.

Ожидается, что в этой новой матрице A’ строки, соответствующие сигналам от одного источника, будут располагаться преимущественно последовательно, что позволит максимально просто группировать их в соответствии с источником соответствующего импульсного шума.

3.2 Применение алгоритма корреляционного анализа

В ходе предварительной обработки экспериментальных данных, было обнаружено, что в миссии 45-46, в которой происходила синхронизация дополнительных седьмого и восьмого каналов БПАС с первым и вторым каналами АР, было обнаружено, что во время длительных сеансов, когда АР и, соответственно, каналы 7 и 8 БПАС были отключены, происходил сбой каналов 6 и 7 БПАС: отключенный канал 7 вместо белого шума выдавал слабый, но «живой» сигнал, а по каналу 6 наблюдалось зашумливание сигнала белым шумом, что заметно сказывалось на высоких частотах (свыше 4 кГц). Во время коротких же сеансов верификации, когда включались оба блока АР и БПАС, все каналы показывали истинные данные.

Из-за этого для единообразия обработки было принято решение проводить корреляционный анализ по первым пяти датчикам БПАС. Датчик М5 находится в той же приборной зоне, что и датчик М6 (таблица 2, рисунки 10 и 11), к тому же в удалении от рабочих мест и кают космонавтов и АСУ, поэтому для координатного разделения источников достаточно лишь одного М5.

Таким образом, программа корреляционного анализа была запущенна на 17×5=85-мерном векторном пространстве. Множитель 5 соответствует количеству датчиков, а множитель 17 - количеству 1/3-октавных полос частот с центральными частотами от 400 Гц включительно до 16 кГц включительно.

На рисунке 22 представлена та же матрица, что и на рисунке 19, но преобразованная посредством одного прохода линейного алгоритма, разработанного в работе [15].

Рис. 22. Упорядоченная линейным алгоритмом корреляционная матрица

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

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

Рис. 23. Упорядоченная рекурсивным алгоритмом корреляционная матрица

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

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

Заключение

В ходе работы была исследована научная информация, полученная научной аппаратурой «Средства определения координат точки пробоя» в ходе космического эксперимента «Пробой» с декабря 2014 года по ноябрь 2015 года.

Получены спектральные характеристики шумящей аппаратуры, действующей на борту РС МКС.

Подтверждено, что НА «СОКП» помимо регистрации и оперативного определения координат точки пробоя может вести мониторинг акустической обстановки на борту МКС.

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

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

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

Список использованных источников

1.      Патент 2387966 РФ. Устройство определения координат места пробоя гермооболочки пилотируемого космического объекта и способ определения координат места пробоя / Авершьев С.П. и др. - Заяв. 16.10.2008; Опубл. 27.04.2010г.

.        Патент 2387965 РФ. Устройство определения координат места пробоя гермооболочки непилотируемого космического объекта и способ определения координат места пробоя / Болотин В.А. и др. - Заяв. 16.10.2008; Опубл. 27.04.2010г.

.        Maeda, N. (1985). A method for reading and checking phase times in autoprocessing system of seismic wave data, Zisin=Jishin 38, 365-379.

.        David P. Schaff, Gotz H. R. Bokelmann, William L. Ellsworth, Eva Zanzerkia, Felix Waldhauser, and Gregory C. Beroza, Optimizing correlation techniques for improved earthquake location, Bull. Seism. Soc. Am., vol. 94, no. 2, pp. 705-721, 2004.

.        Андреева И.Б., Бреховских Л.М. Акустика океана. URL: http://www.akin.ru/spravka/s_ocean.htm (дата обращения 25.02.2016).

.        Shazam: алгоритмы распознавания музыки, сигнатуры, обработка данных. URL: https://habrahabr.ru/company/wunderfund/blog/275043/ (дата обращения 08.02.2016).

.        Как Яндекс распознаёт музыку с микрофона. URL: https://habrahabr.ru/company/yandex/blog/181219/ (дата обращения 08.02.2016).

8.      How does Shazam work. URL: <http://coding-geek.com/how-shazam-works/> (дата обращения 08.02.2016).

.        William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery. 2.6 Singular Value Decomposition //Numerical Recipes in C. - 2nd ed. - Cambridge: Cambridge University Press. -ISBN 0-521-43108-5.

.        Jolliffe I.T. Principal Component Analysis, Series: Springer Series in Statistics, 2nd ed., Springer, NY, 2002, XXIX, 487 p. 28 illus. ISNB 978-0-387-95442-4

.        Жуков А.И. Метод Фурье в вычислительной математике. М. Наука. Физ-матлит. 1992.

.        Половнев А.Л. Программа-анализатор для определения спектральных и корреляционных характеристик случайных виброакустических процессов по данным телеметрических измерений. // Сборник публикаций ЦАГИ. Проблемы аэрокосмической науки и техники. - 2003. - №3.

.        Haijiang Zhang, Clifford Thurber, and Charlotte Rowe, Automatic p-wave arrival detection and picking with multiscale wavelet analysis for single-component recordings, Bull. Seism.Soc. Am., vol. 93, no. 5, pp. 1904-1912, October 2003.

.        Akaike, H. (1973). Information theory and an extension of the maximum likelihood principle, in 2nd Intermational Symposium on Information Theory, Budapest Akademiai Kiado, pp. 267-281.

15.    Зайцев К.И. Алгоритм обработки нестационарных акустических процессов, полученных в космическом эксперименте «Пробой» научной аппаратурой «Средства определения координат точки пробоя». Междисциплинарная курсовая работа. МИЭМ НИУ ВШЭ, 2015.

.        Авершьев С.П., Половнев А.Л. Компьютерная программа с организацией потокового ввода данных с АЦП с одновременным расчетом координат вероятного места пробоя корпуса гермоотсека и записью получаемых данных на диск компьютера в реальном масштабе времени. ЦНИИмаш, фонд алгоритмов и программ. № 851.7553682.4116-01.12, 2007.

Приложение 1. Реализация алгоритма Фурье-анализа на С++

#include "windows.h"

#include <formatio.h>

#include <analysis.h>

#include "stdio.h"

#include "math.h"

#include "ansi_c.h"

double XYZ[8][3];xyz[3];GetVariance(double *wave, DWORD size) {i;mean, s;(size<2) return 0;= 0;(i = 0; i < size; i++) mean+=wave[i];/=size;= 0;(i = 0; i < size; i++) {q;= wave[i]-mean;+=(q*q);

}/=(size-1);s;

}GetWaveFrontAIC(double *wave) {

// AIC-picker from Maeda (1985)*x, *AICn, s1, m1, sN, mN, varx1, varxN, min;lN, l1, N, k, i;= wave;= 4096;= (double*) malloc(N*sizeof(double));(x[0]==x[1]) {= 1;(x[0]==x[i] && i<N-1) i++;(i<N) x[0]=x[i];

}= x[0]*x[0];= x[0];= 1;= 0;(k=1; k<N; k++) sN+=(x[k]*x[k]);= 0;(k=1; k<N; k++) mN+=x[k];= N-1;(k=1; k<N-2; k++) {= l1+1;= s1+x[k]*x[k];= m1+x[k];= (s1-m1*m1/l1)/(l1-1);(varx1<=0) {= GetVariance(x, k+1);(varx1==0) varx1 = 1e-10;

}= lN-1;= sN-x[k]*x[k];= mN-x[k];= (sN-mN*mN/lN)/(lN-1);(varxN<=0) {= GetVariance(&x[k+1], N - k - 1);(varxN==0) varxN = 1e-10;

}[k] = k*log10(varx1)+(N-k-1)*log10(varxN);

}[0]=AICn[1];[N-2]=AICn[N-3];[N-1]=AICn[N-2];= AICn[0];= 0;(i = 0; i < N; i++) if (AICn[i]<min) {= i;= AICn[i];

}(AICn);k;

}dlevel (double * dr, double a, double b) {dret = 0.0;aa = (int)floor(a), bb = (int)ceil(b);(int j=aa; j<=bb; j++) {+=dr[j];(j==bb) dret-=dr[j];

}dret;

}main(int argc, char ** argv) {ch = 6;//8;i, j, k, l;*pf0;(argc>4) {= fopen(argv[2], "w");(pf0, "Time");(l=10; l<27; l++) fprintf(pf0, "\tf%i", (int)(pow(10,(double)l/10-1.4)*1000+0.5));

} else pf0 = fopen(argv[2], "a");FH = OpenFile (argv[3], VAL_READ_ONLY, VAL_OPEN_AS_IS, VAL_BINARY);lip[4096][ch];drre[ch][4096], drim[ch][4096], drf, drg, car6[2049];drre_copy[8][4096], tn[8];[0][0]=-2.667;   XYZ[0][1]=0.850;         XYZ[0][2]=-0.465;[1][0]=-2.637;    XYZ[1][1]=-0.780;         XYZ[1][2]=0.435;[2][0]=-5.785;     XYZ[2][1]=-0.028;         XYZ[2][2]=-0.720;[3][0]=-5.737;    XYZ[3][1]=0.050;         XYZ[3][2]=0.720;[4][0]=-8.926;     XYZ[4][1]=-0.330;         XYZ[4][2]=-0.720;[5][0]=-8.098;    XYZ[5][1]=0.463;         XYZ[5][2]=0.720;[6][0]=0.0;         XYZ[6][1]=0.490;         XYZ[6][2]=-0.571;[7][0]=-0.045;    XYZ[7][1]=-0.757;         XYZ[7][2]=0.156;time[100];(FH, -100*sizeof(SYSTEMTIME), 2);(FH, (char*) time, 100*sizeof(SYSTEMTIME));ntime[100];[0]=0;(i=1; i<100; i++){= time[i].wHour*3600000+time[i].wMinute*60000+time[i].wSecond*1000+time[i].wMilliseconds

time[i-1].wHour*3600000-time[i-1].wMinute*60000-time[i-1].wSecond*1000-time[i-1].wMilliseconds;[i]=0;(j < 153 && j > 0) ntime[i-1]=1;

}(FH, 107, 0);(i=0; i<100; i++) {(FH, (char*)lip, ch*4096*sizeof(int));(ntime[i]==0) {(j=0; j<4096; j++) for(k=0; k<ch; k++) {[k][j] = (double) lip[j][k];_copy[k][j] = drre[k][j];[k][j] = 0.0;(j<2049)car6[j]=0;

}= 0.0;= 0.0;(k=0; k<ch; k++) {(drre[k], drim[k], 4096);(j=0; j<2049; j++) {[k][j]=sqrt(drre[k][j]*drre[k][j]+drim[k][j]*drim[k][j])/2048/(1+(j==0||j==2049));[j]+=drre[k][j];

}+=log10(dlevel(drre[k],30,2049));

}((time[i].wMonth>2 && time[i].wMonth<11)||(time[i].wMonth==2 && time[i].wDay>26)) {[0][0]=-3.845;  XYZ[0][1]=0.998;         XYZ[0][2]=-0.465;

} else {[0][0]=-2.667;    XYZ[0][1]=0.850;         XYZ[0][2]=-1.055;

}(pf0, "\n%d.%d.%d %d:%d:%d_%d", time[i].wYear, time[i].wMonth, time[i].wDay,[i].wHour, time[i].wMinute, time[i].wSecond, time[i].wMilliseconds);(k=0;k<ch;k++) tn[k]=(double)GetWaveFrontAIC(drre_copy[k])/50000.0;(l=10; l<27; l++) fprintf(pf0, "\t%f",(dlevel(car6, pow(10,(double)l/10-1.45)/50*4096, pow(10,(double)l/10-1.35)/50*4096)));

}

}(FH);(pf0);

}

Приложение 2. Реализация алгоритма корреляционного анализа на языке С#

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using System.Threading;

using System.Globalization;System.IO;System.Drawing;proboj

{Program

{static double[] xmass;static double[,] darr;static byte nullsumm(byte[,] dcorr, int i)

{nsum = 0;(int j = 0; j < i; j++) for (int k = i; k < dcorr.GetLength(0); k++) if (nsum < dcorr[j, k]) nsum = dcorr[j, k];nsum;

}static double[,] chij(double[,] dcorr, int i, int j)

{k;t;(k = 0; k < dcorr.GetLength(0); k++)

{= dcorr[k, i];[k, i] = dcorr[k, j];[k, j] = t;

}(k = 0; k < dcorr.GetLength(0); k++)

{= dcorr[i, k];[i, k] = dcorr[j, k];[j, k] = t;

}dcorr;

}static byte[,] chij(byte[,] dcorr, int i, int j)

{k;t;(k = 0; k < dcorr.GetLength(0); k++)

{= dcorr[k, i];[k, i] = dcorr[k, j];[k, j] = t;

}(k = 0; k < dcorr.GetLength(0); k++)

{= dcorr[i, k];[i, k] = dcorr[j, k];[j, k] = t;

}dcorr;

}static double[,] chi(double[,] arr, int i, int j)

{k;t;(k = 0; k < arr.GetLength(1); k++)

{= arr[i, k];[i, k] = arr[j, k];[j, k] = t;

}arr;

}static string getargs(string sorce, int count)

{c = sorce.IndexOf("\t");(c < 0)

{(count == 0) return sorce;return "";

}

{(count == 0) return sorce.Substring(0, c);return getargs(sorce.Substring(c + 1), count - 1);

}

}static void corrmxy(double[,] dcorr, int CMin, int CMax, double diff, string[] time)

{(CMax == 1) return;(diff <= 0) return;i, j, k;sss;tx, ty;nxy;[,] dcorr2 = new byte[CMax, CMax];(i = 0; i < CMax; i++)

{.WriteLine("CorrT " + diff.ToString() + " " + i.ToString());(j = i; j < CMax; j++)

{= 0;(k = 0; k < CMax; k++) nxy += (byte) (Math.Floor(dcorr[i + CMin, k + CMin] + diff) * Math.Floor(dcorr[j + CMin, k + CMin] + diff));[i, j] = nxy;(i != j) dcorr2[j, i] = dcorr2[i, j];

}

}(i = 1; i < CMax - 1; i++)

{.WriteLine("RangC " + diff.ToString() + " " + i.ToString());(j = i + 1; j < CMax; j++)

{= 0;= 0;(k = 0; k < i; k++)

{+= dcorr2[k, i];+= dcorr2[k, j];

}(ty > tx)

{= time[i + CMin];[i + CMin] = time[j + CMin];[j + CMin] = sss;= chij(dcorr, i + CMin, j + CMin);= chij(dcorr2, i, j);= chi(darr, i + CMin, j + CMin);

}

}

}cxl = 0;(i = 1; i < CMax; i++)

{.WriteLine("Diag " + diff.ToString() + " " + i.ToString());(nullsumm(dcorr2, i) == 0)

{(dcorr, cxl + CMin, i - cxl, diff - .01, time);= i;(xmass[i + CMin] == 0) xmass[i + CMin] = diff;

}

}= null;(dcorr, cxl + CMin, CMax - cxl, diff - .01, time);

}void Main(string[] args)

{int fmax = 17*5;double diff = 0.2;[] lines = System.IO.File.ReadAllLines(@"E:\VL_Temp\output-2g6.txt");CMax = lines.Length - 1;i, j, k;lx, ly, lxy;[] time = new string[CMax];];[] loglevel = new double[fmax];= new double[CMax, fmax];[,] dcorr = new double[CMax, CMax];= new double[CMax];(j = 0; j < fmax; j++) loglevel[j] = 0;(i = 0; i < CMax; i++)

{.WriteLine("Read " + i.ToString());[i] = getargs(lines[i + 1], 0);(j = 0; j < fmax; j++)

{[i, j] = Convert.ToDouble(getargs(lines[i + 1], j + 1));[j] += Math.Log10(darr[i, j]);

}

}(j = 0; j < fmax; j++) loglevel[j] /= CMax;(i = 0; i < CMax; i++)

{.WriteLine("Corr " + i.ToString());(j = i; j < CMax; j++)

{= ly = lxy = 0;(k = 0; k < fmax; k++)

{tci = Math.Log10(darr[i, k]);tcj = Math.Log10(darr[j, k]);+= (tci - loglevel[k]) * (tcj - loglevel[k]);+= (tci - loglevel[k]) * (tci - loglevel[k]);+= (tcj - loglevel[k]) * (tcj - loglevel[k]);

}[i, j] = (lxy / Math.Sqrt(lx * ly) + 1) / 2;(i != j) dcorr[j, i] = dcorr[i, j];

}

}bmp = new Bitmap(CMax, CMax);sw = System.IO.File.AppendText(@"E:\VL_Temp\Lex++" + diff.ToString() + ".txt");(i = 0; i < CMax; i++) sw.Write("\t" + time[i]);(i = 0; i < CMax; i++)

{.WriteLine("Out0 " + i.ToString());.Write("\n" + time[i]);(j = 0; j < CMax; j++)

{.Write("\t" + (dcorr[i, j]).ToString());.SetPixel(i, j, Color.FromArgb((int)(255 - 255 * dcorr[i, j] * dcorr[i, j]),

(int)(255 - 255 * dcorr[i, j] * dcorr[i, j]), (int)(255 - 255 * dcorr[i, j] * dcorr[i, j])));

}

}.Close();.Save(@"E:\VL_Temp\Bex++" + diff.ToString() + ".bmp");(dcorr, 0, CMax, diff, time);= new Bitmap(CMax, CMax);= System.IO.File.AppendText(@"E:\VL_Temp\Hex++" + diff.ToString() + ".txt");(i = 0; i < CMax; i++) sw.Write("\t" + time[i]);(i = 0; i < CMax; i++)

{.WriteLine("Out " + i.ToString());.Write("\n" + time[i]);(j = 0; j < CMax; j++)

{.Write("\t" + Math.Floor(100 * dcorr[i, j]).ToString());.SetPixel(i, j, Color.FromArgb((int)(255 - 255 * dcorr[i, j] * dcorr[i, j]),

(int)(255 - 255 * dcorr[i, j] * dcorr[i, j]), (int)(255 - 255 * dcorr[i, j] * dcorr[i, j])));

}(j = 0; j < fmax; j++) sw.Write("\t" + darr[i, j].ToString());.Write("\t" + xmass[i].ToString());

}.Close();.Save(@"E:\VL_Temp\Jex++" + diff.ToString() + ".bmp");

}

}

}