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

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

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

Ожидается, что в этой новой матрице 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г.

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

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

4. 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.

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

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

7. Как Яндекс распознаёт музыку с микрофона. 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).

9. 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.

10. 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

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

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

13. 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.

14. 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.

16. Авершьев С.П., Половнев А.Л. Компьютерная программа с организацией потокового ввода данных с АЦП с одновременным расчетом координат вероятного места пробоя корпуса гермоотсека и записью получаемых данных на диск компьютера в реальном масштабе времени. ЦНИИмаш, фонд алгоритмов и программ. № 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];

double xyz[3];

double GetVariance(double *wave, DWORD size) {

int i;

double mean, s;

if (size<2) return 0;

mean = 0;

for (i = 0; i < size; i++) mean+=wave[i];

mean/=size;

s = 0;

for (i = 0; i < size; i++) {

double q;

q = wave[i]-mean;

s+=(q*q);

}

s/=(size-1);

return s;

}

int GetWaveFrontAIC(double *wave) {

// AIC-picker from Maeda (1985)

double *x, *AICn, s1, m1, sN, mN, varx1, varxN, min;

int lN, l1, N, k, i;

x = wave;

N = 4096;

AICn = (double*) malloc(N*sizeof(double));

if (x[0]==x[1]) {

i = 1;

while (x[0]==x[i] && i<N-1) i++;

if (i<N) x[0]=x[i];

}

s1 = x[0]*x[0];

m1 = x[0];

l1 = 1;

sN = 0;

for (k=1; k<N; k++) sN+=(x[k]*x[k]);

mN = 0;

for (k=1; k<N; k++) mN+=x[k];

lN = N-1;

for (k=1; k<N-2; k++) {

l1 = l1+1;

s1 = s1+x[k]*x[k];

m1 = m1+x[k];

varx1 = (s1-m1*m1/l1)/(l1-1);

if (varx1<=0) {

varx1 = GetVariance(x, k+1);

if (varx1==0) varx1 = 1e-10;

}

lN = lN-1;

sN = sN-x[k]*x[k];

mN = mN-x[k];

varxN = (sN-mN*mN/lN)/(lN-1);

if (varxN<=0) {

varxN = GetVariance(&x[k+1], N - k - 1);

if (varxN==0) varxN = 1e-10;

}

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

}

AICn[0]=AICn[1];

AICn[N-2]=AICn[N-3];

AICn[N-1]=AICn[N-2];

min = AICn[0];

k = 0;

for (i = 0; i < N; i++) if (AICn[i]<min) {

k = i;

min = AICn[i];

}

free(AICn);

return k;

}

double dlevel (double * dr, double a, double b) {

double dret = 0.0;

int aa = (int)floor(a), bb = (int)ceil(b);

for (int j=aa; j<=bb; j++) {

dret+=dr[j];

if(j==bb) dret-=dr[j];

}

return dret;

}

void main(int argc, char ** argv) {

int ch = 6;//8;

int i, j, k, l;

FILE *pf0;

if(argc>4) {

pf0 = fopen(argv[2], "w");

fprintf(pf0, "Time");

for(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");

int FH = OpenFile (argv[3], VAL_READ_ONLY, VAL_OPEN_AS_IS, VAL_BINARY);

int lip[4096][ch];

double drre[ch][4096], drim[ch][4096], drf, drg, car6[2049];

double drre_copy[8][4096], tn[8];

XYZ[0][0]=-2.667; XYZ[0][1]=0.850; XYZ[0][2]=-0.465;

XYZ[1][0]=-2.637; XYZ[1][1]=-0.780; XYZ[1][2]=0.435;

XYZ[2][0]=-5.785; XYZ[2][1]=-0.028; XYZ[2][2]=-0.720;

XYZ[3][0]=-5.737; XYZ[3][1]=0.050; XYZ[3][2]=0.720;

XYZ[4][0]=-8.926; XYZ[4][1]=-0.330; XYZ[4][2]=-0.720;

XYZ[5][0]=-8.098; XYZ[5][1]=0.463; XYZ[5][2]=0.720;

XYZ[6][0]=0.0; XYZ[6][1]=0.490; XYZ[6][2]=-0.571;

XYZ[7][0]=-0.045; XYZ[7][1]=-0.757; XYZ[7][2]=0.156;

SYSTEMTIME time[100];

SetFilePtr (FH, -100*sizeof(SYSTEMTIME), 2);

ReadFile(FH, (char*) time, 100*sizeof(SYSTEMTIME));

int ntime[100];

ntime[0]=0;

for(i=1; i<100; i++){

j = 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;

ntime[i]=0;

if (j < 153 && j > 0) ntime[i-1]=1;

}

SetFilePtr (FH, 107, 0);

for(i=0; i<100; i++) {

ReadFile (FH, (char*)lip, ch*4096*sizeof(int));

if (ntime[i]==0) {

for(j=0; j<4096; j++) for(k=0; k<ch; k++) {

drre[k][j] = (double) lip[j][k];

drre_copy[k][j] = drre[k][j];

drim[k][j] = 0.0;

if(j<2049)car6[j]=0;

}

drf = 0.0;

drg = 0.0;

for(k=0; k<ch; k++) {

FFT(drre[k], drim[k], 4096);

for(j=0; j<2049; j++) {

drre[k][j]=sqrt(drre[k][j]*drre[k][j]+drim[k][j]*drim[k][j])/2048/(1+(j==0||j==2049));

car6[j]+=drre[k][j];

}

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

}

if ((time[i].wMonth>2 && time[i].wMonth<11)||(time[i].wMonth==2 && time[i].wDay>26)) {

XYZ[0][0]=-3.845; XYZ[0][1]=0.998; XYZ[0][2]=-0.465;

} else {

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

}

fprintf(pf0, "\n%d.%d.%d %d:%d:%d_%d", time[i].wYear, time[i].wMonth, time[i].wDay,

time[i].wHour, time[i].wMinute, time[i].wSecond, time[i].wMilliseconds);

for(k=0;k<ch;k++) tn[k]=(double)GetWaveFrontAIC(drre_copy[k])/50000.0;

for(l=10; l<27; l++) fprintf(pf0, "\t%f",

log10(dlevel(car6, pow(10,(double)l/10-1.45)/50*4096, pow(10,(double)l/10-1.35)/50*4096)));

}

}

CloseFile(FH);

fclose(pf0);

}



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

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using System.Threading;

using System.Globalization;

using System.IO;

using System.Drawing;

namespace proboj

{

class Program

{

public static double[] xmass;

public static double[,] darr;

public static byte nullsumm(byte[,] dcorr, int i)

{

byte nsum = 0;

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

return nsum;

}

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

{

int k;

double t;

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

{

t = dcorr[k, i];

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

dcorr[k, j] = t;

}

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

{

t = dcorr[i, k];

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

dcorr[j, k] = t;

}

return dcorr;

}

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

{

int k;

byte t;

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

{

t = dcorr[k, i];

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

dcorr[k, j] = t;

}

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

{

t = dcorr[i, k];

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

dcorr[j, k] = t;

}

return dcorr;

}

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

{

int k;

double t;

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

{

t = arr[i, k];

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

arr[j, k] = t;

}

return arr;

}

public static string getargs(string sorce, int count)

{

int c = sorce.IndexOf("\t");

if (c < 0)

{

if (count == 0) return sorce;

else return "";

}

else

{

if (count == 0) return sorce.Substring(0, c);

else return getargs(sorce.Substring(c + 1), count - 1);

}

}

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

{

if (CMax == 1) return;

if (diff <= 0) return;

int i, j, k;

string sss;

double tx, ty;

byte nxy;

byte[,] dcorr2 = new byte[CMax, CMax];

for (i = 0; i < CMax; i++)

{

Console.WriteLine("CorrT " + diff.ToString() + " " + i.ToString());

for (j = i; j < CMax; j++)

{

nxy = 0;

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

dcorr2[i, j] = nxy;

if (i != j) dcorr2[j, i] = dcorr2[i, j];

}

}

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

{

Console.WriteLine("RangC " + diff.ToString() + " " + i.ToString());

for (j = i + 1; j < CMax; j++)

{

tx = 0;

ty = 0;

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

{

tx += dcorr2[k, i];

ty += dcorr2[k, j];

}

if (ty > tx)

{

sss = time[i + CMin];

time[i + CMin] = time[j + CMin];

time[j + CMin] = sss;

dcorr = chij(dcorr, i + CMin, j + CMin);

dcorr2 = chij(dcorr2, i, j);

darr = chi(darr, i + CMin, j + CMin);

}

}

}

int cxl = 0;

for (i = 1; i < CMax; i++)

{

Console.WriteLine("Diag " + diff.ToString() + " " + i.ToString());

if (nullsumm(dcorr2, i) == 0)

{

corrmxy(dcorr, cxl + CMin, i - cxl, diff - .01, time);

cxl = i;

if (xmass[i + CMin] == 0) xmass[i + CMin] = diff;

}

}

dcorr2 = null;

corrmxy(dcorr, cxl + CMin, CMax - cxl, diff - .01, time);

}

static void Main(string[] args)

{

const int fmax = 17*5;

const double diff = 0.2;

string[] lines = System.IO.File.ReadAllLines(@"E:\VL_Temp\output-2g6.txt");

int CMax = lines.Length - 1;

int i, j, k;

double lx, ly, lxy;

string[] time = new string[CMax];];

double[] loglevel = new double[fmax];

darr = new double[CMax, fmax];

double[,] dcorr = new double[CMax, CMax];

xmass = new double[CMax];

for (j = 0; j < fmax; j++) loglevel[j] = 0;

for (i = 0; i < CMax; i++)

{

Console.WriteLine("Read " + i.ToString());

time[i] = getargs(lines[i + 1], 0);

for (j = 0; j < fmax; j++)

{

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

loglevel[j] += Math.Log10(darr[i, j]);

}

}

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

for (i = 0; i < CMax; i++)

{

Console.WriteLine("Corr " + i.ToString());

for (j = i; j < CMax; j++)

{

lx = ly = lxy = 0;

for (k = 0; k < fmax; k++)

{

double tci = Math.Log10(darr[i, k]);

double tcj = Math.Log10(darr[j, k]);

lxy += (tci - loglevel[k]) * (tcj - loglevel[k]);

lx += (tci - loglevel[k]) * (tci - loglevel[k]);

ly += (tcj - loglevel[k]) * (tcj - loglevel[k]);

}

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

if (i != j) dcorr[j, i] = dcorr[i, j];

}

}

Bitmap bmp = new Bitmap(CMax, CMax);

StreamWriter sw = System.IO.File.AppendText(@"E:\VL_Temp\Lex++" + diff.ToString() + ".txt");

for (i = 0; i < CMax; i++) sw.Write("\t" + time[i]);

for (i = 0; i < CMax; i++)

{

Console.WriteLine("Out0 " + i.ToString());

sw.Write("\n" + time[i]);

for (j = 0; j < CMax; j++)

{

sw.Write("\t" + (dcorr[i, j]).ToString());

bmp.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])));

}

}

sw.Close();

bmp.Save(@"E:\VL_Temp\Bex++" + diff.ToString() + ".bmp");

corrmxy(dcorr, 0, CMax, diff, time);

bmp = new Bitmap(CMax, CMax);

sw = System.IO.File.AppendText(@"E:\VL_Temp\Hex++" + diff.ToString() + ".txt");

for (i = 0; i < CMax; i++) sw.Write("\t" + time[i]);

for (i = 0; i < CMax; i++)

{

Console.WriteLine("Out " + i.ToString());

sw.Write("\n" + time[i]);

for (j = 0; j < CMax; j++)

{

sw.Write("\t" + Math.Floor(100 * dcorr[i, j]).ToString());

bmp.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])));

}

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

sw.Write("\t" + xmass[i].ToString());

}

sw.Close();

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

}

}

}

Размещено на Allbest.ru