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

Проведенное моделирование показало, что подавление помех выполняется уровнемне менее 25 дБ. Также было предложено решение проблемы оценки корреляционной матрицы помех,путем увеличения выборки обучающих пакетов. Предложены методы углового сверхразрешения, необходимые для определения угловых координат источников радиоизлучения. Проведено моделирование их эффективности определения координат в программной среде программирования Matlab.Разработка структурной схемы цифровой адаптивной антенной решетки. Произведен анализ технических параметров предъявляемых цифровой аппаратуре для реализации активной фазированной антенной решетки с цифровым формированием диаграммы направленности. Список использованной литературыРатынский М.В. Адаптация и сверхразрешение в антенных решетках. М.: Радио и связь, 2003. 200 с.Джиган В.И. Адаптивная фильтрация сигналов: теория и алгоритмы. М.: Техносфера, 2013. 528 с. Van Trees H. Optimum Array Processing. New York: Wiley-Interscience, 2002. 1400 pp.Воскресенский Д.И., Канащенкова А.И. Активные фазированные антенные решетки. М.: Радиотехника, 2004. 490 с.Сергиенко А.Б. Цифровая обработка сигналов. СПб.: Питер, 2002. 608 с.Грант П.М., Коуэн К.Ф.Н. Адаптивные фильтры. М.: Мир, 1988. 392 с.Van Veen B.D., Buckley K.M. Beamforming: A versatile approach to spatial filtering. IEEE ASSP Magazine, Vol. 5 No. 2 pp. 4 -24.Чистяков В.А. Определение угловых координат радиоизлучения в антенных решетках космических аппаратов. Universum: Технические науки: электрон. научн. журн. 2017. №12(45).Marple S. Lawrence. Digital Spectral Analysis, Englewood Cliffs, NJ: PrenticeHall, 1987. 373 pp.Зотов С.А., Макаров Е.С., Нечаев Ю.Б. Методы сверхразрешения в задачах радиопеленгации. М.: Известия вузов, 2008. 9 с.Ермолаев В.Т., Флаксман А.Г. Методы оценивания параметров источников сигналов и помех, принимаемых антенной решеткой. Изд.: ННГУ, 2007. 100 с.Чистяков В.А., Куприц В.Ю. Сверхразрешение в антенных решетках: матер. конф. Изд.: СФУ, 2016. 310-312 с.Федорков Б.Г., Телец В.А. Микросхемы ЦАП и АЦП. Функционирование, параметры, применение. М.: Энергоатомиздат, 1990. 320 с.Бахтиаров Г.Д., Малинин В.В., Школин В.П. Аналого-цифровые преобразователи. Проектирование электронной аппаратуры на интегральных микросхемах. М.: Советское радио, 1980. 280 с.Каталог элементов фирмы AnalogDevices: http://www.analog.com/media/en/news-marketing-collateral/product-selection-guide/Space_Products_Selection_Guide.pdfКаталог элементов фирмы TexasInstrument: http://www.ti.com/lsds/ti/space-high-reliability/data-converter/analog-to-digital-converters-products.page#p1498=SpaceМоисеев А.В., Новиков И.Г. Информационный сборник. Книга 1.Каталог элементов фирмы Microsemi:http://www.actel.ru/rubric/radiacionno-stojkie-plisКаталог элементов фирмы Xilinx:https://www.xilinx.com/products/silicon-devices/fpga.htmlКомаров А.В. Цифровые сигнальные процессоры. Обнинск, 2003. 141 с.Приложение АЛистинг алгоритма адаптивной обработки сигналовclear all;close all;i=sqrt(-1);N=25; % кол-воэлементовФАРf=8e9; % частотаГГцc = physconst('LightSpeed'); % скоростьсветаlamda=c/f; % длинаволныd=lamda/2; % расстояние между элементами ФАР%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% формирование плоской ФАР %%for n=1:5 for m=1:5 r(m,n)=d*(n-1); % формирование элементов с расстоянием между ними d=lamda/2 endendpl_az=reshape(r,[1,25]); % расположение элементов для азимутальной плоскости (азимут)pl_ug=reshape(r',[1,25]); % расположение элементов для плоскости возвышения (угол места)%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% формирование узкополосных полезных сигналов %% Fd=125e3;k=1:Fd;Az1=5/57.3; % азимут сигнала 1 (10 градусов)Ug1=0/57.3; % угол места сигнала 1 (0 градусов)f1=35e3; % частота полезного сигнала 1A1=1e-5; % мощность сигнала 1 в разахs1=A1*sin(2*pi*k*f1/Fd);Az2=0/57.3; % азимут сигнала 2 (0 градусов)Ug2=-7/57.3; % угол места сигнала 2 (-17 градусов)f2=30e3; % частота полезного сигнала 2A2=1e-5; % мощность сигнала 2 в разахs2=A2*sin(2*pi*k*f2/Fd);%%%%%% формирование набега фаз узкополосных полезных сигналов %%%%%% %%%%%% сигнал 1 %%%%%%az_array_sig1=pl_az*sin(Az1); % набег фазы сигнала 1 по азимуту ug_array_sig1=pl_ug*sin(Ug1); % набег фазы сигнала 1 по углу местаsig1=az_array_sig1+ug_array_sig1; % набег фазы сигнала 1 по всей ФАРarray_sig1=exp(sig1*(-2*i*pi/lamda)); % набег фазы сигнала 1 с учетом расположения элементов по всей ФАРS1=s1.*array_sig1'; % сигнал 1 пришедший на вход ФАР с учетом набега фаз%%%%%% сигнал 2 %%%%%%az_array_sig2=pl_az*sin(Az2); % набег фазы сигнала 2 по азимуту ug_array_sig2=pl_ug*sin(Ug2); % набег фазы сигнала 2 по углу местаsig2=az_array_sig2+ug_array_sig2; % набег фазы сигнала 2 по всей ФАРarray_sig2=exp(sig2*(-2*i*pi/lamda)); % набег фазы сигнала 2 с учетом расположения элементов по всей ФАРS2=s2.*array_sig2'; % сигнал 2 пришедший на вход ФАР с учетом набега фаз%%%%%%%%%%%%%%%%%%%% формирование широкополосных помех %%Az_p1=30/57.3; % азимут помехи 1 (30 градусов)Ug_p1=0/57.3; % угол места помехи 1 (0 градусов)E1=1e-2; % мощность помехи 1 (разы)p1=E1*(randn(1,125000)+i*randn(1,125000));Az_p2=0/57.3; % азимут помехи 2 (0 градусов)Ug_p2=-45/57.3; % угол места помехи 2 (-45 градусов)E2=1e-2; % мощность помехи 2 (разы)p2=E2*(randn(1,125000)+i*randn(1,125000));Az_p3=-60/57.3; % азимут помехи 3 (-60 градусов)Ug_p3=0/57.3; % угол места помехи 3 (0 градусов)E3=1e-2; % мощность шума (разы)p3=E3*(randn(1,125000)+i*randn(1,125000));Az_p4=0/57.3; % азимут помехи 4 (0 градусов)Ug_p4=36/57.3; % угол места помехи 4 (36 градусов)E4=1e-2; % мощность шума (разы)p4=E4*(randn(1,125000)+i*randn(1,125000));%%%%%% формирование набега фаз широкополосных помех %%%%%% %%%%%% помеха 1 %%%%%%az_array_p1=pl_az*sin(Az_p1); % набег фазы помехи 1 по азимуту ug_array_p1=pl_ug*sin(Ug_p1); % набег фазы помехи 1 по углу местаpom1=az_array_p1+ug_array_p1; % набег фазы помехи 1 по всей ФАРarray_p1=exp(pom1*(-2*i*pi/lamda)); % набег фазы помехи 1 с учетом расположения элементов по всей ФАРP1=p1.*array_p1'; % помеха 1 пришедшая на вход ФАР с учетом набега фаз%%%%%% помеха 2 %%%%%%az_array_p2=pl_az*sin(Az_p2); % набег фазы помехи 2 по азимуту ug_array_p2=pl_ug*sin(Ug_p2); % набег фазы помехи 2 по углу местаpom2=az_array_p2+ug_array_p2; % набег фазы помехи 2 по всей ФАРarray_p2=exp(pom2*(-2*i*pi/lamda)); % набег фазы помехи 2 с учетом расположения элементов по всей ФАРP2=p2.*array_p2'; % помеха 2 пришедшая на вход ФАР с учетом набега фаз%%%%%% помеха 3 %%%%%%az_array_p3=pl_az*sin(Az_p3); % набег фазы помехи 3 по азимуту ug_array_p3=pl_ug*sin(Ug_p3); % набег фазы помехи 3 по углу местаpom3=az_array_p3+ug_array_p3; % набег фазы помехи 3 по всей ФАРarray_p3=exp(pom3*(-2*i*pi/lamda)); % набег фазы помехи 3 с учетом расположения элементов по всей ФАРP3=p3.*array_p3'; % помеха 3 пришедшая на вход ФАР с учетом набега фаз%%%%%% помеха 4 %%%%%%az_array_p4=pl_az*sin(Az_p4); % набег фазы помехи 4 по азимуту ug_array_p4=pl_ug*sin(Ug_p4); % набег фазы помехи 4 по углу местаpom4=az_array_p4+ug_array_p4; % набег фазы помехи 4 по всей ФАРarray_p4=exp(pom4*(-2*i*pi/lamda)); % набег фазы помехи 4 с учетом расположения элементов по всей ФАРP4=p4.*array_p4'; % помеха 4 пришедшая на вход ФАР с учетом набега фаз%%%%%% сумма всех сигналов на выходе ФАР %%%%%SIG=(S1+S2+P1+P2+P3+P4); % сумма на 25 элементахSIG_array=sum(SIG); % сумма на выходе ФАР, проссумированных на 25 элементах%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Алгоритм адаптивной фильтрации %%%%%P=P1+P2+P3+P4; % сумма помех поступающих на вход антенной решеткиKor_matr=(P*P')+diag(ones(1,25)); % оценка корреляционной матрицы помехObr_Kor_matr=inv(Kor_matr); % оценка обратной корреляционной матрицы помех%%%%% формирование управляющего вектора диаграммой направленности %%%%%%Az=0/57.3;Ug=0/57.3;az_array_A=pl_az*sin(Az); ug_array_A=pl_ug*sin(Ug); A=az_array_A+ug_array_A; % набег фазы сигнала 1 по всей ФАР array_A=exp(i*A*(-2*pi/lamda)); % набег фазы сигнала 1 с учетом расположения элементов по всей ФАР%%%%%% формирование веткора весовых коэффициентов для адаптивной фильтрации помехW_opt=Obr_Kor_matr*array_A';%%%%% Построение диаграммы направленности ФАР после адаптации %%%%%Array=phased.URA([5 5],[lamda/2, lamda/2]);%%%%% ДН с адаптацией %%%%%plotResponse(Array,f,c,'RespCut','Az','Format','Line','Weights',W_opt); % ДН азимутальной плоскостиfigure;plotResponse(Array,f,c,'RespCut','El','Format','Line','Weights',conj(W_opt)); % ДН плоскости возвышения figure;plotResponse(Array,f,c,'RespCut','Az','Format','Polar','Weights',conj(W_opt)); %ДН азимутальной плоскости (полярная)figure;plotResponse(Array,f,c,'RespCut','El','Format','Polar','Weights',conj(W_opt)); %ДН плоскости возвышения (полярная)figure;plotResponse(Array,f,c,'RespCut','3d','Format','Line','Weights',conj(W_opt)); %ДН 3D figure;plotResponse(Array,f,c,'RespCut','3d','Format','Polar','Weights',conj(W_opt)); %ДН 3D (полярная)%{%%%%% скрипт для вывода ДН с адаптацией и без на одном графике %%%%%% для использования удалить фигурные скобки { и }plotResponse(Array,f,c,'RespCut','Az','Format','Line','Weights',W_opt); % ДН азимутальной плоскостиhold on;plotResponse(Array,f,c,'RespCut','Az','Format','Line');% ДН азимутальной плоскости без адаптацииfigure;plotResponse(Array,f,c,'RespCut','El','Format','Line','Weights',conj(W_opt)); % ДН плоскости возвышения hold on;plotResponse(Array,f,c,'RespCut','El','Format','Line'); % ДН плоскости возвышения без адаптации%}%%%%% Спектральный анализатор для проверки адаптации %%%%%SIG1=SIG'*W_opt; % входной сигнал (адитивная смесь полезных сигналов и помех) после адаптацииsignalAnalyzer(real(SIG_array),'SampleRate',125e3); % входной сигнал до адаптации signalAnalyzer(real(SIG1),'SampleRate',125e3); % входной сигнал после адаптации%signalAnalyzer(real(s1+s2)','SampleRate',125e3); % полезные сигналы %signalAnalyzer(real(sum(P1+P2+P3+P4))','SampleRate',125e3); % помехиПриложение БЛистинг программы углового сверхразрешенияclear all;close all;i=sqrt(-1);%%%%%% параметры плоской линейной фазированной антенной решетки (ФАР) %%%%%%N=16; % кол-во элементов ФАРf=8e9; % частота ГГцc = physconst('LightSpeed'); % скорость светаlamda=c/f; % длина волны d=lamda/2; % расстояние между элементами ФАРfor n=1:N r(n)=d*(n-1); % формирование элементов с расстоянием между ними d=lamda/2end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% формирование сигналов % узкополосный сигнал 1Fd=125e3;k=1:Fd;Az1=15/57.3; % азимут сигнала 1 (10 градусов)f1=35e3; % частота полезного сигнала 1A1=1e4; % мощность сигнала 1 в разахs1=A1*sin(2*pi*k*f1/Fd);az_array_s1=r*sin(Az1);s1_az=exp(az_array_s1*(-2*i*pi/lamda));S1=s1.*s1_az';% узкополосныйсигнал 2Fd=125e3;k=1:Fd;Az2=-27/57.3; % азимут сигнала 2 (-27 градусов)f2=15e3; % частота полезного сигнала 2A2=1e4; % мощность сигнала 2 в разахs2=A2*sin(2*pi*k*f2/Fd);az_array_s2=r*sin(Az2);s2_az=exp(az_array_s2*(-2*i*pi/lamda));S2=s2.*s2_az';% широкополосный сигнал 3Az_s3=30/57.3; % азимут сигнала 3 (30 градусов)A3=1e4; % мощность сигнала 3 (разы)s3=A3*(randn(1,125000)+i*randn(1,125000));az_array_s3=r*sin(Az_s3); % набегфазысигнала 3s3_az=exp(az_array_s3*(-2*i*pi/lamda)); % набег фазы сигнала 3 с учетом расположения элементов по всей ФАРS3=s3.*s3_az'; % помеха 3 пришедшая на вход ФАР с учетом набега фаз% широкополосный сигнал 4Az_s4=-70/57.3; % азимут сигнала 4 (30 градусов)A4=1e4; % мощность сигнала 4 (разы)s4=A4*(randn(1,125000)+i*randn(1,125000));az_array_s4=r*sin(Az_s4); % набегфазысигнала 4s4_az=exp(az_array_s4*(-2*i*pi/lamda)); % набег фазы сигнала 4 с учетом расположения элементов по всей ФАРS4=s4.*s4_az'; % помеха 3 пришедшая на вход ФАР с учетом набега фазS=S1+S2+S3+S4; % сумма всех сигналов на входе ФАРKor_matr=(S*S')+diag(ones(1,N))/(A1+A2+A3+A4); % оценкакорреляционнойматрицыпомех%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Алгоритмсверхразрешения (поискасигналов) MUSIC[Sbv,Sbc]=eig(Kor_matr); % % Гистограмма для определения числа сигналовHH=diag(Sbc)';CC=1:N;bar(CC,HH);Z=4; % число столбцов (это и есть сигналы)En=[Sbv(:,1:(N-Z))];% Z-число собственных чисел сигнального пространства, которые нужно откинуть (число сигналов)% построение пеленгационного рельефаNQ=100000;%Шаг сканированияm = 1:1:NQ; gr=-90+180./NQ*m;for m2=1:NQ tetta = -pi/2+ pi/NQ*m2; for m3=1:N V(m3,1)=exp(i*pi*(m3-1)*sin(tetta)); end V1=V'; QS(m2)=(1/(V1*En*En'*V));end max11=max(abs(QS));diag11=10*log10(abs(QS)/max11);figure;plot(gr,diag11);title('Алгоритм MUSIC');xlabel('Угловые координаты источника радиоизлучения, гр.');ylabel('Пеленгационный рельеф, дБ.');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Алгоритм сверхразрешения (поиска сигналов) Кейпона%вычисление обратной корреляционной матрицы сигналов OKM=inv((S*S')+diag(ones(1,N)));for m2=1:NQ tetta = -pi/2+ pi/NQ*m2; for m3=1:N V(m3,1)=exp(i*pi*(m3-1)*sin(tetta)); end V1=V'; Q(m2)=1/(V1*OKM*V); endmax1=max(abs(Q));diag1=10*log10(abs(Q)/max1);figure;plot(gr,diag1);title('АлгоритмКейпона');xlabel('Угловые координаты источника радиоизлучения, гр.');ylabel('Пеленгационный рельеф, дБ.');