数字波束形成递推最小二乘算法研究+源代码.docx

数字波束形成递推最小二乘算法研究摘要：数字波束形成(DBF)是一门新兴的阵列天线技术，自适应波束形成(ADBF)针对复杂信号环境对波束形成的一个很好的补充。自适应技术的核心思想是优化理论，常用的最佳准则有最小均方误差(MMSE)，最大信噪比(MaxSNR)，线性约束最小方差(LCMv)，最大似然(ML)，最小二乘法(LS)等；本问主要对递推最小二乘算法(RLS)进行了包括了流程设计、性能分析和仿真评价。关键字：DBF、ADBF、RLS、一 数字波束形成简介数字波束形成 (DBF)1是随着数字信号处理方法的发展而建立起来的一门新 技术，这种技术不仅能充分保留阵列天线上收集的信息，而且能利用复杂的数字 信号信号处理方法对信息进行处理，因此它具有一般雷达不具有的超低分辨率和 低副瓣的性能，波束扫描，自适应波束控制等。由于采用了先进的数字信号处理 方法和自适应技术，DBF雷达不仅性能优越，而且非常灵活，被认为是下代雷达 系统中广泛使用的一项新技术。相对于一般模拟波束形成技术，数字波束形成技术具有非常多的优点，比如 超低分辨率，超低副瓣，抗电子干扰，扫描迅速，多目标处理，高性能并行数字 处理等。自适应波束形成(ADBF)针对复杂信号环境对波束形成的一个很好的补充。 它的核心思想是采用各种优化算法和自适应算法，根据阵列对不同环境的响应， 产生不同的加权系数，己取得最好的信号特征。自适应波束形成可以实现自适应 的空一时处理，进行各种非线性处等，波束灵活，抗干扰强，具有较好的自校正和 自适应能力。阵列信号处理是现代信号处理2的一个分支，其本质是利用空间分散排列的传 感器阵列和多通道接收机来获取信号的时域和空域等多维信息，以达到检测信号 和提取其参数的目的。阵列信号处理的主要内容可分为波束形成技术，空间谱估计等方面技术，他们都是基于对信号进行空间采样的的数据进行处理，因此这些技术是相互渗透和相互关联的。二 数字波束形成原理和典型波束形成准则2.1 数字波束形成原理假设接收天线为N阵元均匀直线阵3，目标的来波方向为，第k个波束指向为k，相邻阵元间距为d，信号波长为，且各阵元都是各向同性的，对K阵元的加权分别为w1k，w2k,wNK，信号是窄带信号。简单地讲，数字多波束形成器就是一个乘加器，如图2-1所示。表示第i阵的接收信号经过A/D变换和数字正交后的复数字信号。加权系数，其中表示提供的阵内相位补偿值，表示降低天线副瓣所需的幅度加权。进行相位和幅度补偿后，并对各阵列的输出信号相加，即可得： (2-1)对式(2-1)进行求模，就可以得到第K个波束的天线方向图函数。天线的方向图不是唯一的，根据要求，相同的数据可以用不同方法加权(改变权因子)，以便形成不同形式的波束和任意多的波束。通过数字波束合成后输出的信号(如信号确实投射在该波束内)可以大幅度的提高信号的SNR。最后得到的复输出信号再直接送入后续处理单元。图2-1 自适应数字波束形成原理图表示成矩阵形式如下：(其中wk为权矢量，Xk为数据矢量，T表示转置，它们都为列矢量。2.2 自适应数字波束形成的原理DBF技术是当代雷达的热点技术之一，对于阵列信号处理来说，它往往和自 适应技术相互联系。因为面对非常复杂的环境时候，恶化的工作环境将会使不具 备自我适应，自我调整能力的DBF雷达的性能急剧下降。因此，为了提高DBF系 统的抗干扰能力，必须求助于自适应数字波束形成(ADBF)技术，自适应阵列是实现ADBF的基础。自适应阵列的结构框图如图2-1示。从图中可以看出自适应阵列是由按一定空间排列的多个阵单元构成，是一种在实际环境下自行控制其方向图的天线系统。它能实时地对外界未知的干扰环境作出反应，在干扰的到达方向形成零点或降低此方向的副瓣电平，这样便可以保证 接收所需信号与干扰噪声的信噪比有一个最佳值。自适应阵列正是利用这种空间 特性，改善了阵列输出SNR，抑制了强干扰。自适应技术的核心思想是优化理论4，我们熟悉的优化理论有拉格朗日求极值函数的变分法，最陡下降法，最小二乘法等。 目前常用的最佳准则有最小均方误差(MMSE)，最大信噪比(MaxSNR)，线性约束最小方差(LCMv)，最大似然(ML)，最小二乘法(LS)等，下面简单介绍几种较常使用的准则。2.3 最小均方误差准则(MMSE)最小均方误差准则就是使阵列输出(K代表K时刻)与参考信号d(k)的均方误差最小，均方误差为： (2-2)其中r=E(d(k)Xk),R=E(XXH),一般地将R称为互相关矩阵5。将式（2-4）对于权向量求梯度，得到梯度算子： 令梯度算子为零，得到最小均方误差准则下的最佳权向量：2.4 最大信噪比原则（MaxSNR）最大信噪比原则是基于期望信号的功率与噪声功率之比的最大的准则，假设期望信号为S，且Rt=E（ssH），RW=E(uuH),其中u表示噪声，则有：这时，输出的信噪比SNR为：经计算，使得输出信噪比最大的最佳权向量是对应于矩阵的最大特征值的特征值的特征向量，得到的最佳权系数满足：2.5 线性约束最小方差准则（LCMV）在已知期待信号的来波方向和参考信号的条件下，最小方差准则是通过最小化阵列输出的噪声方差，来取得对信号S的较好的增益。经权重后的波束形成器的输出为: 为保证波束形成对信号S增益，必须对波束形成器的权向量加以限制，使其在信号S的方向产生一定的增益，即：其中，为期望信号的方向矢量，则最佳权重可以表示为：三 递推最小二乘算法（RLS）3.1 自适应算法简介自适应算法所采用的最优准则有最小均方误差(LMS)准则6，最小二乘(LS)准则、最大信嗓比准则和统计检测准则等，其中最小均方误差(LMS)准则和递推最小二乘法(LS)准则是目前最为流行的自适应算法准则。LMS算法和RLS算法由于采用的最优准则不同，因此这两种算法在性能，复杂度等方面均有许多差别。对于某个自适应算法，衡量其算法性能的指标主要为以下几个：（1）算法的收敛速度。用算法达到稳定状态 (即与最优值的接近程度达到一定值) 需要迭代的次数表示。（2）算法的跟踪性能。它指信道发生变化时，算法自适应跟踪信道的能力。（3）算法的稳健性。只当输入病态的情况下算法能否正常工作。（4）算法计算的复杂程度。自适应波束形成算法可有多种分类方式，一般的分类方法根据算法是否需要显式的训练序列，可分为非盲算法和盲算法两类，另外，根据权值作用于阵元或波束，还可分为阵元空间和波束空间两类。盲算法一般利用调制信号本身固有的、与具体承载的信息比特无关的一些特征，如恒模 、子空间、有限符号集和高阶统计量等。与非盲算法相比，它无需发送已知符号，故可节省系统频谱资源。其缺点是一般收敛速度较慢，且存在相位模糊问题。非盲算法相对盲算法而言，通常误差较小，收敛速度也较快，但需要占用一定的系统资源用于传送参考信号、训练序列等。即使非盲算法收敛速度快，但仍然跟不上快衰落变化的速率要求。对于通信而言，算法的收敛速度是很重要的一个性质。所以本文主要介绍非盲算法中的RLS算法。非盲算法主要包括最小均方误差算法LMS、采样矩阵求逆法SMI和迭代最小二乘算法RLS。LMS 算法由于其收敛速度很慢，因此在信号环境变化很快时其性能不是很好，且算法性能对阵列信号协方差矩阵的特征值散布度很敏感，当散布度很大时，很难收敛。 SMI算法收敛速度比LMS快，但是算法复杂了不少，而且在硬件实现中的有限字长效应会给求逆运算带来数值上的不稳定。RLS算法是基于使每一快拍的阵列输出平方和最小的准则，即最小二乘（LS）准则。它利用了从算法初始化后得到的所有阵列数据信息，用递推方法来完成矩阵的求逆运算，因而收敛速度快，对特征值的散布度不敏感，且能实现收敛速度与计算复杂性之间的折衷。一般在大信噪比的情况下，RLS 比LMS的收敛速度快一个数量级。3.2递推最小二乘算法(RLS)RLS算法是SMI算法的迭代过程，该算法的性能准则是使参考信号和阵列输出信号之间的加权均方误差最小。与LMS不同的是，RLS算法的性能量度准则是使信号的加权平方误差和最小化。RLS算法在权值迭代调整过程中用一个与信号相关的附加矩阵代替LMS算法中的步长，从而提高了收敛速度。RLS算法的权向量迭代式为： （3-17） （3-18）其中表示参考信号，为迭代过程中的过渡矩阵，其中称为遗忘因子()。引入遗忘因子的作用是让离n时刻较近的误差有较大的权重，距离较远的拥有降低的权重，确保以前观测到的数据被渐渐“遗忘”，从而使滤波器工作在一个平稳状态下。下面是算法推导：设：则：令：则：其中 为增益向量，又： 所以式中：先验误差3.3 总结RLS算法的步骤1、初始化：w(0)=0, R(0)=I,2、更新：对于n=1、2计算：希望相关矩阵初始值R(0)在R(n)中占很小的比重，因此设R(0)=I。一般取0.001。对RLS算法的总体评价：1、RLS算法对非平稳信号的适应性好。2、RLS算法收敛速度快，估计精度高稳定性好。3、遗忘因子越大，越不易遗忘，效果越好。4、RLS算法计算复杂度高，不利于实时性处理。四 递推最小二乘法(RLS)仿真结果及分析仿真实验使用的是的Matlab软件，MATLAB是由美国Math Works公司推出的用于数值计算和图形处理的科学计算系统环境。MATLAB是英文MATrix LABoratory(矩阵实验室)的缩写，它集中了日常数学处理中的各种功能，包括高效的数值计算、矩阵运算、信号处理和图形生成功能。MATLAB是一个功能强大的系统，它集数值计算、矩阵运算、信号处理和图形生成功能为一体，另外，MATLAB还具有很强的功能扩展能力，可以配置各种各样的工具箱，以完成一些特定的任务，同时，用户还可以根据自己的工作任务，开发自己的工具箱。在MATLAB环境下，可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。4.1 递推最小二乘法(RLS)算法性能仿真分析图4.1 图4.2图4.1和图4.2所示为RLS在不同的遗忘因子取值下的收敛速度曲线。图4.1为=0.4时的收敛速度，图4.2所示的是=0.98时的收敛速度。从曲线中可以看出遗忘因子的取值越大，收敛速度越慢。图 4.3图 4.4图4.3和图4.4所示为RLS在不同的遗忘因子取值下误差的估计值。图4.3为=0.4时的误差估计值，图4.4所示的是=0.98时的误差估计值。图片显示遗忘因子越大，误差越大，稳定性越差。如果遗忘因子很小，则信号对期望信号的依赖性就会很大，所以，输出信号就很接近期望信号，这就导致了如果期望信号是错误的，那么输出信号也是错误的。所以综合考虑，遗忘因子应该取在0.95和0995之间。4.2 递推最小二乘法(RLS)应用于数字波束形成的仿真分析按照给定要求，天线阵元数16个，辐射源5个，在遗传因子为0.95，预设有用信号角度为30度，干扰信号角度分别为-80，0，60，80.预子图误差的对比111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111首先是单次运算均方误差和100次均方误差的仿真对比图图 4.4由图4.4可见，在单次和100次运算下的均方误差相差不大。图 4.5由图4.5可见，该算法生成的方向图满足要求，在覆盖180度的方位面内，产生了一个波束，波束指向为30度，为有用信号角度；同时对干扰方向-80，0 ，60，80都有一定程度的置零；另外对杂波信号进行了压制，主瓣波束变窄。图 4.6图4.6表示的是输入输出信号的对比图，可以见到通过该算法有效的还原了有用信号，并且对噪声和干扰信号都有了较好的抑制。五 总结自适应数字波束形成算法在阵列天线技术中起着至关重要的做用，本文通过先介绍了数字波束形成简介及其原理，后比较了几个常用的数字波束形成准则，然后引入了递推最小二乘算法(RLS)，对其进行算法描述和推导，最后对算法性能和应用于算法的性能进行了仿真分析。六 参考文献 1吴伟陵.移动通信中的关键技术M.北京邮电大学出版社，2000:256-268.2石新智,王高峰.修正MUSCI算法对非线性阵列适用性的讨论.电子学报J,2004,32(1) :15-18.3龚耀寰编著.自适应滤波时域自适应滤波和智能天线M.北京电子工业出版社，2003:217-307.4刘君,廖桂生.信源数目过估计和欠估计下MUSIC算法分析J.2004,26(2)：123-137.5 Harry LVan Trees.Optimum array processing.Publication by John Wiley and sons LtdM.2002:710-916.6S Bellofiore,J F Balanis.Smart-Antenna Systems for mobile Comm Networks Park Beamforming andNetwork Throughput. IEEE antennas and Propagation MagazineJ.2002,44(4):106-114. 附录1算法性能测试代码：M=15;%均衡滤波器阶数为2*M+1Lb=10;%信道b长度为L+1%Lb=2;%信道c长度为L+1%hb=0.4070.8150.407;%离散时间信道chb=0.04-0.050.07-0.210.500.720.360.000.210.030.07;%离散时间信道bHb=zeros(2*M+1,2*M+Lb+1);fork=1:2*M+1;%信道b的信道矩阵 Hb(k,k:1:k+Lb)=hb;end%产生伯努利序列和加性白噪声，构建均衡滤波器的输入数据矩阵sigma=1e-3;%加性白高斯噪声的方差N=2000;%迭代次数s=randsrc(2*M+Lb+N,1);%伯努利序列vn=sqrt(sigma)*randsrc(2*M+Lb+N,1);S=zeros(2*M+Lb+1,N);%发射信号矩阵SV=zeros(2*M+1,N);%加性白高斯噪声矩阵Vfork=1:NS(:,k)=s(2*M+Lb+k:-1:k);V(:,k)=vn(2*M+k:-1:k);endUb=Hb*S+V;%均衡滤波器输入数据矩阵Ub%RLS迭代算法dn=S(M+Lb+1,:);%期望信号lambda=0.990;%RLS遗忘因子delta=0.004;%RLS调整参数wb_RLS=zeros(2*M+1,N+1);wb_RLS(M+1,1)=1;%权向量初始值epsilon=zeros(N,1);%先验估计误差P1=eye(2*M+1)/delta;%相关矩阵逆的初始值fork=1:N%RLS算法迭代过程 PIn=P1*Ub(:,k);deno=lambda+Ub(:,k)*PInkn=PIn/deno;epsilon(k)=dn(k)-wb_RLS(:,k)*Ub(:,k);wb_RLS(:,k+1)=wb_RLS(:,k)+kn*conj(epsilon(k);P1=P1/lambda-kn*Ub(:,k)*P1/lambda;endMSEB_RLS=abs(epsilon).2;%单次实验均方误差MSEB_RLSn=1:2000;plot(n,MSEB_RLS)title(RLS算法学习曲线)xlabel(迭代次数n)ylabel(MSEB_RLS)axis(020001e-31e+2);附录2波束形成算法代码：Clc,clear all,close all;M=16;%阵元数%d=bochang/2;%阵元间距为信号波长的一半N=1000;%快拍数theta_signal=30;%有用信号的方向theta_interfere1=0;%干扰信号1的方向theta_interfere2=80;%干扰信号2的方向theta_interfere3=-80;%干扰信号3的方向theta_interfere4=60;%干扰信号4的方向snr=10;%所有信号的信噪比均为此值linear_snr=10(snr/10);%线性信噪比值t=1:N; f0=2e7;%参考信号的频率s0=5*cos(2*pi*f0*t/(N*f0);%目标信号for i=1:M fai0(i,1)=exp(j*(i-1)*pi*sin(pi*theta_signal/180);%有用信号的相位 fai1(i,1)=exp(j*(i-1)*pi*sin(pi*theta_interfere1/180);%干扰信号1的相位 fai2(i,1)=exp(j*(i-1)*pi*sin(pi*theta_interfere2/180);%干扰信号2的相位 fai3(i,1)=exp(j*(i-1)*pi*sin(pi*theta_interfere3/180);%干扰信号3的相位 fai4(i,1)=exp(j*(i-1)*pi*sin(pi*theta_interfere4/180);%干扰信号2的相位endpower=0;for t=1:N power=power+s0(t)2;endpower_signal=power/N;%有用信号功率power_noise=power_signal/linear_snr;%噪声功率power_interfere1=power_noise*linear_snr;%干扰信号1的功率power_interfere2=power_noise*linear_snr;%干扰信号2的功率power_interfere3=power_noise*linear_snr;%干扰信号2的功率power_interfere4=power_noise*linear_snr;%干扰信号2的功率s1=normrnd(0,sqrt(power_interfere1/2),1,N)+j*normrnd(0,sqrt(power_interfere1/2),1,N);%干扰信号1的随机包络s2=normrnd(0,sqrt(power_interfere2/2),1,N)+j*normrnd(0,sqrt(power_interfere2/2),1,N);%干扰信号2的随机包络s3=normrnd(0,sqrt(power_interfere3/2),1,N)+j*normrnd(0,sqrt(power_interfere3/2),1,N);%干扰信号3的随机包络s4=normrnd(0,sqrt(power_interfere4/2),1,N)+j*normrnd(0,sqrt(power_interfere4/2),1,N);%干扰信号4的随机包络X0=; noise=tco_wgn(M,N,0.2,2);for t=1:N x=s0(t)*fai0+s1(t)*fai1+s2(t)*fai2+s3(t)*fai3+s4(t)*fai4+noise(:,t); %阵列对信号的完整响应,x为M*1 X0=X0;x;%X0为N*Mend%*RLS算法*%w=zeros(M,1);%权矢量初始化lamda=1;%遗忘因子lamdadelta=10;c=delta*eye(M,M);%接收信号样本自相关矩阵的逆矩阵初始化pingjunzhi=zeros(1,N); mse1=;for n=1:N miu=X0(n,:)*