基于矩阵分解的多输入多输出雷达解相关算法

基于波束空间MUSIC的MIMO雷达波达方向估计作者：盛志超盛骥松来源：《现代电子技术》2011年第17期摘要：在进行波达方向估计时，阵元空间MUSIC方法的计算量通常都比较大。

为了解决此问题，采用了波束空间MUSIC的方法，它的计算量较阵元空间MUSIC方法有所下降，将它运用于多输入多输出雷达波达方向的估计问题。

计算机仿真实验表明，虽然协方差矩阵特征分解的计算量下降了，但是波束空间MUSIC的性能依然良好。

关键词：MIMO雷达；波达方向估计；波束空间； MUSIC中图分类号：TN958-34 文献标识码：A文章编号：1004-373X(2011)17-0018-03MIMO Radar DOA Estimation Based on Beam-space MUSICSHENG Zhi-chao1, SHENG Ji-song2(1. Jiangsu University of Science and Technology, Zhenjiang 212003, China; 2. The 723 Institute of CSIC, Yangzhou 225001, China)Abstract： The array-space MUSIC requires large calculated amount usually during DOA estimation. In order to resolve this problem, the algorithm based on beam-space MUSIC is adopted. The calculated amount of this algorithm is lower than array-space MUSIC and is used in MIMO radar DOA estimation. Simulation results show that the performance of beam-space MUSIC is good all the same though the calculated amount of eigen-decomposition of covariance matrix is reduced.Keywords： MIMO radar; DOA estimation; beam-space; MUSIC0 引言近年来，MIMO雷达成为国内外雷达领域中研究的热点。

mimo矩阵分解
MIMO (multiple-input multiple-output) 矩阵分解是一种针对MIMO 系统中的问题进行矩阵分解的方法，它将 MIMO 系统中的通信信道矩阵分解为两个相互独立的矩阵，分别表示发送和接收端之间的信号处理过程。

这种矩阵分解可用于信道估计、信号检测、波束形成等应用中。

MIMO 矩阵分解的基本思想是将 MIMO 信道矩阵分解为两个矩阵：一个发射矩阵和一个接收矩阵。

发射矩阵表示发送端的信号处理，接收矩阵表示接收端的信号处理。

通过这种矩阵分解，可以将 MIMO 信道模型简化为多个单输入单输出信道模型，从而简化信号处理过程。

同时，该矩阵分解技术可以大幅减少信道估计和信号检测中需要的计算复杂度。

MIMO 矩阵分解有不同的方法，包括奇异值分解 (SVD)、QR 分解、Cholesky 分解等。

其中，SVD 是最常用的方法之一，可以将信道矩阵分解为三个矩阵，包括一个左奇异矩阵、一个右奇异矩阵和一个奇异值矩阵。

这种分解方法可以保证最小化信道矩阵的条件数，并提高通信系统的性能。

总之，MIMO 矩阵分解是一种用于 MIMO 信道中信号处理和通信系统设计的重要技术，可以简化信号处理过程，提高通信系统的性能和效率。

多输入多输出系统(MIMO ,Multiple input multiple output)最早是控制系统中提出的一个概念，它表示一个系统有多个输入和多个输出。

而MIMO技术早期用于干扰无线信号，后来则用于移动通信和固定宽度的无线领域。

如果将移动通信系统的传输信道看成一个系统，则发射信号可看成移动信道(系统)的输入信号，而接收信号则可看成移动信道(系统)的输出信号。

在通信中，由多径引起的衰落通常被认为是有害因素，不过对于MIMO系统而言，多径可引起的衰落以作为一个有利因素并加以利用。

MIMO 技术以其可以有效利用多径引起的衰落来成倍地提高业务传输速率，并引发了通信的一次革命。

基于通信系统中的MIMO技术的使用情况，近几年国外学者提出了MIMO雷达的概念。

1.MIMO雷达信号处理发展历史1.1 国外研究现状国外最早在MIMO雷达信号处理领域开始开创性的工作者有New Jersey Institute of Technology的Eran Fishler、Alex Haimovich等人，他们研究的工作主要集中在MIMO雷达的信号建模，从模型中获取我们感兴趣的参数的算法研究（如散射点的散射系数，散射点距雷达的距离等），并从雷达对目标检测性能等方面说明它相对于普通的相控阵雷达所具有的优越性，明确指出了MIMO雷达将是未来雷达发展的一个趋势。

几乎在同一时期，MIT Lincoln Laboratory的K. W. Forsythe等人的研究工作也在同步进行，他们的研究工作则主要集中在MIMO雷达的性能优越性的理论证明。

同时该实验室的Frank C. Robey也作了大量的实验，通过大量实验证明MIMO雷达相比传统的雷达有许多优点。

目前国外研究MIMO雷达的著名机构有美国的佛罗尼达大学(University of Florida)，MIT Lincoln Laboratory、新泽西理工学院（New Jersey Institute of Technology）等。

互耦条件下MIMO雷达非圆目标稳健角度估计方法王咸鹏;国月皓;黄梦醒;沈重;曹春杰;冯文龙【摘要】提出了一种在互耦条件下基于酉张量分解的多输入多输出(MIMO)雷达非圆目标稳健的角度估计算法.所提算法首先在张量域利用互耦系数矩阵的带状对称Toeplitz结构来消除未知互耦的影响,然后通过构造一个特殊的增广张量捕获非圆信号的非圆特性与其固有的多维结构特性,并利用增广张量的centro-Hermitian 特性通过酉变换转化为实值张量,最后利用高阶奇异值分解(HOSVD)获得信号子空间,结合实值子空间技术获得目标的离开方向(DoD)和到达方向(DoA)估计.由于同时利用信号的非圆结构与多维结构特性,所提算法具有比现有的子空间算法更准确的角度估计性能,同时所提算法只需要实值运算,具有较低的运算复杂度.仿真结果表明,所提算法具有有效性与优越性.【期刊名称】《通信学报》【年(卷),期】2019(040)007【总页数】7页(P144-150)【关键词】双基地MIMO雷达;角度估计;非圆信号;互耦;高阶奇异值分解【作者】王咸鹏;国月皓;黄梦醒;沈重;曹春杰;冯文龙【作者单位】海南大学南海海洋资源利用国家重点实验室,海南海口 570228;海南大学信息与科学技术学院,海南海口 570228;海南大学南海海洋资源利用国家重点实验室,海南海口 570228;海南大学信息与科学技术学院,海南海口 570228;海南大学南海海洋资源利用国家重点实验室,海南海口 570228;海南大学信息与科学技术学院,海南海口 570228;海南大学南海海洋资源利用国家重点实验室,海南海口570228;海南大学信息与科学技术学院,海南海口 570228;海南大学南海海洋资源利用国家重点实验室,海南海口 570228;海南大学信息与科学技术学院,海南海口570228;海南大学南海海洋资源利用国家重点实验室,海南海口 570228;海南大学信息与科学技术学院,海南海口 570228【正文语种】中文【中图分类】TN9581 引言多输入多输出（MIMO, multiple input multiple output）雷达概念一经提出，就立即引起了雷达研究领域学者们的广泛关注[1-2]。

第４６卷　第４期２０２４年４月系统工程与电子技术ＳｙｓｔｅｍｓＥｎｇｉｎｅｅｒｉｎｇａｎｄＥｌｅｃｔｒｏｎｉｃｓＶｏｌ．４６　Ｎｏ．４Ａｐｒｉｌ２０２４文章编号：１００１ ５０６Ｘ（２０２４）０４ １４５６ １０　网址：ｗｗｗ．ｓｙｓ ｅｌｅ．ｃｏｍ收稿日期：２０２２１２１３；修回日期：２０２３０６２５；网络优先出版日期：２０２３０８１０。

网络优先出版地址：ｈｔｔｐｓ：∥ｋｎｓ．ｃｎｋｉ．ｎｅｔ／ｋｃｍｓ／ｄｅｔａｉｌ／１１．２４２２．ＴＮ．２０２３０８１０．１１１２．００２．ｈｔｍｌ基金项目：国家自然科学基金（６１６７１０９５，６１７０２０６５，６１７０１０６７，６１７７１０８５）；信号与信息处理重庆市市级重点实验室建设项目（ＣＳＴＣ２００９ＣＡ２００３）；重庆市自然基金（ｃｓｔｃ２０２１ｊｃｙｊｍｓｘｍＸ０８３６）；重庆市教育委员会科研项目（ＫＪ１６００４２７，ＫＪ１６００４２９）资助课题 通讯作者．引用格式：邹涵，张天骐，马 然，等．基于多特征融合的ＭＩＭＯ ＯＦＤＭ系统单混信号调制识别算法［Ｊ］．系统工程与电子技术，２０２４，４６（４）：１４５６ １４６５．犚犲犳犲狉犲狀犮犲犳狅狉犿犪狋：ＺＯＵＨ，ＺＨＡＮＧＴＱ，ＭＡＫＲ，ｅｔａｌ．Ｓｉｎｇｌｅ ｍｉｘｅｄｓｉｇｎａｌｍｏｄｕｌａｔｉｏｎａｎｄｒｅｃｏｇｎｉｔｉｏｎａｌｇｏｒｉｔｈｍｆｏｒＭＩＭＯ ＯＦＤＭｓｙｓｔｅｍｂａｓｅｄｏｎｍｕｌｔｉ ｆｅａｔｕｒｅｆｕｓｉｏｎ［Ｊ］．ＳｙｓｔｅｍｓＥｎｇｉｎｅｅｒｉｎｇａｎｄＥｌｅｃｔｒｏｎｉｃｓ，２０２４，４６（４）：１４５６ １４６５．基于多特征融合的犕犐犕犗 犗犉犇犕系统单混信号调制识别算法邹　涵 ，张天骐，马 然，杨宗方（重庆邮电大学通信与信息工程学院，重庆４０００６５） 摘　要：为解决非协作通信中多输入多输出正交频分复用（ｍｕｌｔｉｐｌｅ ｉｎｐｕｔｍｕｌｔｉｐｌｅ ｏｕｔｐｕｔｏｒｔｈｏｇｏｎａｌｆｒｅ ｑｕｅｎｃｙｄｉｖｉｓｉｏｎｍｕｌｔｉｐｌｅｘｉｎｇ，ＭＩＭＯＯＦＤＭ）系统的单混信号调制识别问题，提出一种基于多特征融合和决策融合的调制识别方法。

ｄｏｉ：１０．３９６９／ｊ．ｉｓｓｎ．１００３－３１１４．２０２２．０２．０１７引用格式：施育鑫，鲁信金，孙艺夫，等．ＭＩＭＯ通信模型下的相关矩阵组低复杂度设计［Ｊ］．无线电通信技术，２０２２，４８（２）：３２７－３３５．［ＳＨＩＹｕｘｉｎ，ＬＵＸｉｎｊｉｎ，ＳＵＮＹｉｆｕ，ｅｔａｌ．Ｌｏｗ⁃ｃｏｍｐｌｅｘｉｔｙＤｅｓｉｇｎｏｆＣｏｒｒｅｌａｔｉｏｎＭａｔｒｉｘＧｒｏｕｐｕｎｄｅｒＭＩＭＯＣｏｍｍｕｎｉｃａｔｉｏｎＭｏｄｅｌ［Ｊ］．ＲａｄｉｏＣｏｍｍｕｎｉｃａｔｉｏｎｓＴｅｃｈｎｏｌｏｇｙ，２０２２，４８（２）：３２７－３３５．］ＭＩＭＯ通信模型下的相关矩阵组低复杂度设计施育鑫１，鲁信金２，孙艺夫２，雷㊀菁２，李玉生１（１．国防科技大学第六十三研究所，江苏南京２１００００；２．国防科技大学电子科学学院，湖南长沙４１００００）摘㊀要：矩阵组常用于无线通信中的数据表示㊂在多输入多输出（Ｍｕｌｔｉｐｌｅ⁃ＩｎｐｕｔＭｕｌｔｉｐｌｅ⁃Ｏｕｔｐｕｔ，ＭＩＭＯ）通信模型中，基站利用信道数据设计适应信道的最小均方误差（ＭｉｎｉｍｕｍＭｅａｎＳｑｕａｒｅＥｒｒｏｒ，ＭＭＳＥ）均衡器接收矩阵组，以低复杂度地处理来自多用户的上行数据㊂首先分析了矩阵数据的关联性，通过时谱图确定矩阵组在时间维所具有的强相关性；其次采用插值算法进行低复杂度的矩阵估计，并提出最大插值比搜索算法计算各类插值算法的性能及其复杂度；接着利用一种改进的Ｓｔｒａｓｓｅｎ矩阵求逆算法来降低ＭＭＳＥ求逆过程的复杂度㊂相比传统的接收矩阵组，显著降低了计算复杂度㊂关键词：ＭＩＭＯ通信模型；奇异值分解；最小均方误差；计算复杂度中图分类号：ＴＮ９２９．５㊀㊀㊀文献标志码：Ａ㊀㊀㊀开放科学（资源服务）标识码（ＯＳＩＤ）：文章编号：１００３－３１１４（２０２２）０２－０３２７－０９Ｌｏｗ⁃ｃｏｍｐｌｅｘｉｔｙＤｅｓｉｇｎｏｆＣｏｒｒｅｌａｔｉｏｎＭａｔｒｉｘＧｒｏｕｐｕｎｄｅｒＭＩＭＯＣｏｍｍｕｎｉｃａｔｉｏｎＭｏｄｅｌＳＨＩＹｕｘｉｎ１，ＬＵＸｉｎｊｉｎ２，ＳＵＮＹｉｆｕ２，ＬＥＩＪｉｎｇ２，ＬＩＹｕｓｈｅｎｇ１（１．６３ｒｄＲｅｓｅａｒｃｈＩｎｓｔｉｔｕｔｅ，ＮａｔｉｏｎａｌＵｎｉｖｅｒｓｉｔｙｏｆＤｅｆｅｎｓｅＴｅｃｈｎｏｌｏｇｙ，Ｎａｎｊｉｎｇ２１００００，Ｃｈｉｎａ；２．ＳｃｈｏｏｌｏｆＥｌｅｃｔｒｏｎｉｃＳｃｉｅｎｃｅ，ＮａｔｉｏｎａｌＵｎｉｖｅｒｓｉｔｙｏｆＤｅｆｅｎｓｅＴｅｃｈｎｏｌｏｇｙ，Ｃｈａｎｇｓｈａ４１００００，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ：Ｍａｔｒｉｘｇｒｏｕｐｓａｒｅｏｆｔｅｎｕｓｅｄｔｏｒｅｐｒｅｓｅｎｔｄａｔａｉｎｗｉｒｅｌｅｓｓｃｏｍｍｕｎｉｃａｔｉｏｎ．ＩｎｔｈｅＭＩＭＯ（Ｍｕｌｔｉｐｌｅ⁃ＩｎｐｕｔＭｕｌｔｉｐｌｅ⁃Ｏｕｔｐｕｔ）ｃｏｍｍｕｎｉｃａｔｉｏｎｍｏｄｅｌ，ｔｈｅｂａｓｅｓｔａｔｉｏｎｏｆｔｅｎｕｓｅｓｔｈｅｃｈａｎｎｅｌｄａｔａｔｏｄｅｓｉｇｎａｃｈａｎｎｅｌ⁃ａｄａｐｔｅｄＭＭＳＥ（ＭｉｎｉｍｕｍＭｅａｎＳｑｕａｒｅＥｒｒｏｒ）ｅｑｕａｌｉｚｅｒｒｅｃｅｉｖｉｎｇｍａｔｒｉｘ，ｉｎｏｒｄｅｒｔｏｐｒｏｃｅｓｓｔｈｅｕｐｌｉｎｋｄａｔａｆｒｏｍｍｕｌｔｉｐｌｅｕｓｅｒｓｗｉｔｈｌｏｗｃｏｍｐｌｅｘｉｔｙ．Ｆｉｒｓｔ，ｔｈｅｃｏｒｒｅｌａｔｉｏｎｏｆｍａｔｒｉｘｄａｔａｉｓａｎａｌｙｚｅｄ，ａｎｄｔｈｅｓｔｒｏｎｇｃｏｒｒｅｌａｔｉｏｎｏｆｔｈｅｍａｔｒｉｘｇｒｏｕｐｉｎｔｈｅｔｉｍｅｄｉｍｅｎｓｉｏｎｉｓｄｅｔｅｒｍｉｎｅｄｔｈｒｏｕｇｈｔｉｍｅ⁃ｓｐｅｃｔｒｏｇｒａｍ．Ｓｅｃｏｎｄｌｙ，ａｎｉｎｔｅｒｐｏｌａｔｉｏｎａｌｇｏｒｉｔｈｍｉｓｕｓｅｄｆｏｒｌｏｗ⁃ｃｏｍｐｌｅｘｉｔｙｍａｔｒｉｘｅｓｔｉｍａｔｉｏｎ，ａｎｄａｍａｘｉｍｕｍｉｎｔｅｒｐｏｌａｔｉｏｎｒａｔｉｏｓｅａｒｃｈａｌｇｏｒｉｔｈｍｉｓｐｒｏｐｏｓｅｄｔｏｃａｌｃｕｌａｔｅｔｈｅｐｅｒｆｏｒｍａｎｃｅａｎｄｃｏｍｐｌｅｘｉｔｙｏｆｖａｒｉｏｕｓｉｎｔｅｒｐｏｌａｔｉｏｎａｌｇｏｒｉｔｈｍｓ．ＴｈｅｎａｎｉｍｐｒｏｖｅｄＳｔｒａｓｓｅｎｍａｔｒｉｘｉｎｖｅｒｓｉｏｎａｌ⁃ｇｏｒｉｔｈｍｉｓｕｓｅｄｔｏｒｅｄｕｃｅｔｈｅｃｏｍｐｌｅｘｉｔｙｏｆｔｈｅＭＭＳＥｉｎｖｅｒｓｉｏｎｐｒｏｃｅｓｓ．Ｃｏｍｐａｒｅｄｗｉｔｈｔｈｅｔｒａｄｉｔｉｏｎａｌｒｅｃｅｉｖｉｎｇｍａｔｒｉｘ，ｔｈｅｃｏｍｐｕｔａ⁃ｔｉｏｎａｌｃｏｍｐｌｅｘｉｔｙｉｓｓｉｇｎｉｆｉｃａｎｔｌｙｒｅｄｕｃｅｄ．Ｋｅｙｗｏｒｄｓ：ＭＩＭＯｃｏｍｍｕｎｉｃａｔｉｏｎｍｏｄｅｌ；ｓｉｎｇｕｌａｒｖａｌｕｅｄｅｃｏｍｐｏｓｉｔｉｏｎ（ＳＶＤ）；ＭＭＳＥ；ｃｏｍｐｕｔａｔｉｏｎａｌｃｏｍｐｌｅｘｉｔｙ收稿日期：２０２１－１２－１６论文来源：基于２０２１年 华为杯 第十八届中国研究生数学建模竞赛Ａ题建模改写，参赛团队获得竞赛一等奖（获奖率约１．１％）及华为专项奖（共１０项）㊂０㊀引言矩阵常常被用于表示无线通信㊁图像视频处理㊁计算机视觉㊁相控阵雷达的数据表示㊂随着用户需求的不断增加，数据规模㊁通信阵列的持续扩大，矩阵的大小和维度也随之快速增加，这给矩阵的数据存储㊁算法计算带来了很大的困难㊂矩阵的关联性是指矩阵数据在某些维度上的相关特性，例如视频中时间相邻的帧具有很强的矩阵关联性㊂因此，充分挖掘矩阵间关联性，以实现低复杂度的计算具有十分重要的价值和意义㊂１㊀问题描述对于所给定复数矩阵Ｈ＝Ｈｊ，ｋ{}，Ｈｊ，ｋɪＭˑＮ，ｊ＝１，２， ，Ｊ；ｋ＝１，２， ，Ｋ㊂其中，矩阵之间以及同一矩阵的元素之间有一定的相关性，包括：相同ｊ下标㊁不同ｋ下标的矩阵间存在一定的关联，即Ｈｊ，１，Ｈｊ，２，Ｈｊ，３， ，Ｈｊ，Ｋ{}间存在关联；且矩阵的各个元素间ｈ（ｊ，ｋ）ｍ，ｎ{}，ｍ＝１，２， ，Ｍ；ｎ＝１，２， ，Ｎ也存在关联，矩阵Ｈｊ，ｋɪＭˑＮ，ｊ＝１，２， ，Ｊ；ｋ＝１，２， ，Ｋ可表示为：Ｈｊ，ｋ＝ｈ（ｊ，ｋ）１，１ｈ（ｊ，ｋ）１，２ｈ（ｊ，ｋ）１，３ｈ（ｊ，ｋ）１，Ｎｈ（ｊ，ｋ）２，１ｈ（ｊ，ｋ）２，２ｈ（ｊ，ｋ）２，３ ｈ（ｊ，ｋ）２，Ｎ︙︙︙︙ｈ（ｊ，ｋ）Ｍ，１ｈ（ｊ，ｋ）Ｍ，２ｈ（ｊ，ｋ）Ｍ，３ ｈ（ｊ，ｋ）Ｍ，Ｎéëêêêêêùûúúúúú㊂（１）此外定义矩阵组Ｈ＝Ｈｊ，ｋ{}的一组数学运算，其中间结果Ｖ＝Ｖｊ，ｋ{}由式（２）给出：Ｖｊ，ｋ＝ｓｖｄ（Ｈｊ，ｋ）Ｈｊ，ｋ＝Ｕｊ，ｋＳｊ，ｋＶ Ｈｊ，ｋＶｊ，ｋ＝Ｖ Ｈｊ，ｋ（：，１：Ｌ）ìîíïïïï，（２）式中，ｊ＝１，２， ，Ｊ；ｋ＝１，２， ，Ｋ，ｓｖｄ（㊃）为矩阵奇异值分解（ＳｉｎｇｕｌａｒＶａｌｕｅＤｅｃｏｍｐｏｓｉｔｉｏｎ，ＳＶＤ）中求解右奇异向量的过程；Ｖｊ，ｋ是由Ｈｊ，ｋ的前Ｌ个右奇异向量构成的矩阵，维度为ＮˑＬ㊂为得到最终输出结果Ｗ＝Ｗｊ，ｋ{}，先将不同ｊ下标㊁相同ｋ下标的Ｖｊ，ｋ进行横向的拼接，得到维度ＮˑＬＪ的Ｖｋ＝［Ｖ１，ｋ， Ｖｊ，ｋ， ＶＪ，ｋ］，然后根据式（３）获取Ｗｋ：Ｗｋ＝Ｖｋ（ＶＨｋＶｋ＋σ２Ｉ）－１，（３）式中，σ２为固定常数；Ｗｋ维度同Ｖｋ；Ｉ为单位矩阵，维度为ＬＪˑＬＪ㊂最后将各Ｗｋ按如式（４）进行拆解：Ｗｋ＝［Ｗ１，ｋ， Ｗｊ，ｋ， ＷＪ，ｋ］，（４）式中，Ｗｊ，ｋ为Ｗｋ中顺序排列的子矩阵，维度为ＮˑＬ㊂为了降低计算和储存的复杂度，分析相关矩阵组的关联性，通过建模对输出结果进行估计，建模过程可用式（５）表示：Ｗ＾＝ｆ（Ｈ），（５）式中，Ｗ＾即为对输出结果Ｗ的建模估计㊂该建模过程可拆分为如式（６）的两个步骤㊂Ｖ＾＝ｆ１（Ｈ）Ｗ＾＝ｆ２（Ｖ＾）{，（６）式中，ｆ１（㊃）表示从输入矩阵组Ｈ到中间结果Ｖ的建模过程，Ｖ＾表示中间结果Ｖ的建模估计，ｆ２（㊃）表示从中间结果Ｖ到最终结果Ｗ的建模过程，Ｗ＾表示最终结果Ｗ的建模估计㊂定义Ｗ的建模估计精度为：ρｌ，ｊ，ｋ（Ｗ）＝Ｗ＾Ｈｌ，ｊ，ｋＷｌ，ｊ，ｋ２Ｗ＾ｌ，ｊ，ｋ２ Ｗｌ，ｊ，ｋ ２，ｌ＝１，２， ，Ｌ，（７）式中， ㊃ ２表示矢量的欧几里得范数（即２范数，对于列矢量ａ， ａ ２＝ａＨａ），Ｗｌ，ｊ，ｋ表示Ｗｊ，ｋ的第ｌ列㊂上式中，Ｗ＾Ｈｌ，ｊ，ｋＷｌ，ｊ，ｋ为复数标量，此处取其欧几里得范数即获取其模值㊂为描述方便，额外定义Ｗ的最低建模精度为ρｍｉｎ（Ｗ）：ρｍｉｎ（Ｗ）ｍｉｎｌɪ｛１，２， ，Ｌ｝ｊɪ｛１，２， ，Ｊ｝ｋɪ｛１，２， ，Ｋ｝ρｌ，ｊ，ｋ（Ｗ），（８）式中，ｍｉｎｌ，ｊ，ｋ（㊃）表示在ｌ，ｊ，ｋ三个维度上取最小值㊂另外，中间结果Ｖ的建模估计精度ρｌ，ｊ，ｋ（Ｗ）的定义及最低建模精度ρｍｉｎ（Ｗ）的定义与此相同㊂计算复杂度定义为由矩阵组Ｈ计算得到结果矩阵组Ｗ所需要的总计算复杂度㊂复数矩阵运算可拆解为基本的复数运算，而基本的复数运算又可进一步拆解为基本的实数运算㊂例如，复数乘法（ａ＋ｂｊ）（ｃ＋ｄｊ）＝（ａｃ－ｂｄ）＋（ａｄ＋ｂｃ）ｊ的复杂度为４次实数乘法和２次实数加（减）法㊂实数基本运算的复杂度按照表１计算㊂表１㊀实数基本运算的计算复杂度Ｔａｂ．１㊀Ｃｏｍｐｕｔａｔｉｏｎａｌｃｏｍｐｌｅｘｉｔｙｏｆｂａｓｉｃｏｐｅｒａｔｉｏｎｓｏｎｒｅａｌｎｕｍｂｅｒｓ运算类型计算复杂度加（减）法１乘法３倒数２５平方根２５自然指数２５自然对数２５正弦２５余弦２５其他１００２０２１研究生数学建模Ａ题提供的数据集附件（Ｄａｔａ１ Ｄａｔａ６）给出的详细数据，包括输入矩阵组㊁标准中间矩阵组和标准输出矩阵组的数据及其维度，其中Ｍ＝４，Ｎ＝６４，Ｌ＝２，Ｊ＝４，Ｋ＝３８４，σ２＝０．０１，数据为十进制格式㊂根据所给数据Ｄａｔａ１ Ｄａｔａ６中的Ｍ＝４，Ｎ＝６４，Ｊ＝４，可对应通信模型中共有Ｊ＝４个用户，每个用户的发射天线数为Ｍ＝４，基站的接收天线数为Ｎ＝６４㊂Ｌ＝２表示取信道衰落程度最低的２个信道向量，即对矩阵进行压缩㊂Ｋ表示信道探测时隙数㊂σ２表示信道中高斯白噪声的噪声方差㊂基于给定的所有矩阵数据，本文通过分析数据间的关联性，解决相关矩阵组的低复杂度计算问题，即以减少计算复杂度为目标进行模型优化㊂设计相应的近似分析模型Ｗ＾＝ｆ（Ｈ），在满足ρｍｉｎ（Ｖ）ȡρｔｈ＝０．９９的情况下，使根据表格计算的总计算复杂度最低㊂２　通信模型建立建立基于矩阵的多输入多输出（Ｍｕｌｔｉｐｌｅ⁃ＩｎｐｕｔＭｕｌｔｉｐｌｅ⁃Ｏｕｔｐｕｔ，ＭＩＭＯ）通信模型［１］如图１所示，Ｊ个用户发送信息，信号经过信道Ｈ到达基站，基站有Ｎ根接收天线㊂图１㊀矩阵关系的通信模型建立示意图Ｆｉｇ．１㊀Ｓｃｈｅｍａｔｉｃｄｉａｇｒａｍｏｆｅｓｔａｂｌｉｓｈｉｎｇｃｏｍｍｕｎｉｃａｔｉｏｎｍｏｄｅｌｏｆｍａｔｒｉｘｒｅｌａｔｉｏｎｓｈｉｐ如图２所示，对于某个用户，拥有Ｍ个天线，各个天线均可与基站天线进行通信㊂其通信信道矩阵Ｈｊ，ｋ＝Ｕｊ，ｋＳｊ，ｋＶＨｊ，ｋ，通过ＳＶＤ分解，将信道矩阵分解成方向酉矩阵和信道随机衰落矩阵的乘积，其中，Ｓｊ，ｋ为随机矩阵，代表波束随机衰落主信道矩阵，Ｕｊ，ｋ和Ｖ Ｈｊ，ｋ分别是用户和基站特征向量矩阵的相关矩阵㊂由于ＶＨｊ，ｋ与基站和用户位置相关且各个节点的位置相对固定，可以取信道衰落程度最低的Ｌ个信道向量压缩Ｖ Ｈｊ，ｋ矩阵，即Ｖｊ，ｋ＝ＶＨｊ，ｋ（：，１：Ｌ）㊂以上模型建立与式（２）一致㊂图２㊀单个用户和基站的通信示意图Ｆｉｇ．２㊀Ｓｃｈｅｍａｔｉｃｄｉａｇｒａｍｏｆｃｏｍｍｕｎｉｃａｔｉｏｎｂｅｔｗｅｅｎａｓｉｎｇｌｅｕｓｅｒａｎｄａｂａｓｅｓｔａｔｉｏｎ此时，可使用压缩矩阵Ｖｊ，ｋ来表示Ｈ矩阵㊂进一步，引入最小均方误差（ＭｉｎｉｍｕｍＭｅａｎＳｑｕａｒｅＥｒｒｏｒ，ＭＭＳＥ）的概念来解释题设条件㊂如图３所示的信道模型中，信号ｘ经过信道Ｖ＝Ｖｊ，ｋ，由于白噪声ｎ的影响，接收信号ｙ可表示为：ｙ＝Ｖｘ＋ｎ㊂（９）图３㊀信号传输模型（求解ＭＭＳＥ流程）Ｆｉｇ．３㊀Ｓｉｇｎａｌｔｒａｎｓｍｉｓｓｉｏｎｍｏｄｅｌ（ｐｒｏｃｅｓｓｏｆｓｏｌｖｉｎｇＭＭＳＥ）ＭＭＳＥ的目的是找到一个矩阵Ｗ＝Ｗｊ，ｋ{}，使得Ｗｙ更加接近ｘ㊂得到ｘ ＝Ｗｙ与原始发送信号ｘ的差值为：ｅ＝ｘ －ｘ＝Ｗｙ－ｘ㊂（１０）此时的ＭＭＳＥ为：ＭＭＳＥ＝ＥｅＨｅ{}㊂（１１）假设接收到的数据ｙ和误差ｅ是不相关的，即Ｅｅ㊃ｙＨ{}＝０㊂（１２）将式（１０）代入式（１２）可得：Ｅ（Ｗｙ－ｘ）㊃ｙＨ{}＝０㊂（１３）将式（１３）左边进一步展开可得：㊀㊀Ｅ（Ｗｙ－ｘ）㊃ｙＨ{}＝ＥＷｙｙＨ{}－ＥｘｙＨ{}＝ＷＥｙｙＨ{}－ＥｘｙＨ{}㊂（１４）由式（１３）和式（１４）可得：Ｗ＝ＥｘｙＨ{}ＥｙｙＨ{}－１㊂（１５）接下来对ＥｙｙＨ{}和ＥｘｙＨ{}进行处理，首先对于ＥｙｙＨ{}，将其进一步展开：㊀㊀ＥｙｙＨ{}＝Ｅ（Ｖｘ＋ｎ）（Ｖｘ＋ｎ）Ｈ{}＝ＥＶｘｘＨＶＨ＋ＶｘｎＨ＋ｎｘＨＶＨ＋ｎｎＨ{}㊂（１６）此处假设输入信号ｘ和噪声ｎ不相关，则ｎｘＨ与ｎｘＨ值为０，可得：㊀㊀㊀ＥｙｙＨ{}＝ＥＶｘｘＨＶＨ＋ｎｎＨ{}＝ＶＥｘｘＨ{}ＶＨ＋ＥｎｎＨ{}＝Ｖ（Ｐ㊃Ｉ）ＶＨ＋σ２㊃Ｉ，（１７）式中，Ｐ为发送信号ｘ的能量，σ２为噪声ｎ的方差㊂其次对于ＥｘｙＨ{}，展开如下：㊀㊀㊀ＥｘｙＨ{}＝Ｅｘ（Ｖｘ＋ｎ）Ｈ{}＝ＥｘｘＨＶＨ＋ｘｎＨ{}＝ＥｘｘＨＶＨ{}＝ＥｘｘＨ{}ＶＨ＝（Ｐ㊃Ｉ）ＶＨ㊂（１８）得到ＥｙｙＨ{}和ＥｘｙＨ{}后，将其代入式（１５）可得到Ｗ的表达式：㊀㊀Ｗ＝ＥｘｙＨ{}ＥｙｙＨ{}－１＝（Ｐ㊃Ｉ）ＶＨ（Ｖ（Ｐ㊃Ｉ）ＶＨ＋σ２㊃Ｉ）－１＝Ｐ㊃ＶＨ（ＰＶＶＨ＋σ２㊃Ｉ）－１＝ＶＨ（ＶＶＨ＋㊃Ｉ）－１㊂（１９）当发送信号ｘ的能量Ｐ为１时，则可得：Ｗ＝ＶＨ（ＶＶＨ＋σ２㊃Ｉ）－１㊂（２０）综合上述分析，ＭＩＭＯ模型中利用信道数据计算信道的ＭＭＳＥ均衡器接收矩阵的复杂度主要来源于ＳＶＤ分解与式（２０）中的矩阵求逆㊂３㊀相关矩阵组的低复杂度计算由前文可知，Ｈ㊁Ｖ和Ｗ之间的关系可以由图４表示㊂图４㊀Ｈ㊁Ｖ和Ｗ的关系示意图Ｆｉｇ．４㊀ＲｅｌａｔｉｏｎｓｈｉｐｏｆＨ，ＶａｎｄＷ３．１㊀利用相关矩阵组的关联性降低计算复杂度利用相关矩阵组的关联性以降低计算复杂度，其具体分析及操作如下㊂３．１．１㊀矩阵数据的关联性对于信道系数复数矩阵Ｈ＝Ｈｊ，ｋ{}，其中Ｈｊ，ｋɪＭˑＮ，ｊ＝１，２， ，Ｊ；ｋ＝１，２， ，Ｋ㊂因此，该矩阵是一个ＭˑＮˑＪˑＫ维的信道矩阵，其中Ｍ表示接收天线的数量，Ｎ表示发射天线的数量，Ｊ表示用户数量，Ｋ表示时隙个数㊂由前文建立的通信模型，对Ｈ矩阵内在的关联性进行分析㊂首先，Ｈｊ，ｋ内部的关系可表示为不同天线构建出的不同信道之间的相关性㊂在一般高斯白噪声信道条件下，天线阵列之间固定的距离和入射角关系将带来一定的规律，但数据集中未能发现Ｈｊ，ｋ内部可靠的相关特性㊂这可能是由于天线之间的方向性㊁距离之间的差异较大，使得信道在空间上的相关特性不再明显㊂考虑不同ｋ时隙，同一用户ｊ的信道系数情况，即ＭＮ个信道的时间相关性㊂图５给出了ＭＮ个信道在时隙ｋ＝１，２，３情况下的幅度响应和相位响应㊂可以看出，在不同ｋ下的幅度和角度的变化程度不大，这可以理解为在相关时间内，信道的变化很小，这进一步验证通信建模的可行性㊂（ａ）不同信道系数的幅度响应（ｂ）不同信道系数的相位响应图５㊀同ｊ不同ｋ的Ｈ块之间的幅度和相位响应关系Ｆｉｇ．５㊀ＡｍｐｌｉｔｕｄｅａｎｄｐｈａｓｅｒｅｓｐｏｎｓｅｒｅｌａｔｉｏｎｓｈｉｐｂｅｔｗｅｅｎｂｌｏｃｋｓＨｗｉｔｈｔｈｅｓａｍｅｊａｎｄｄｉｆｆｅｒｅｎｔｋ基于上述两层分析，得出Ｈｊ，ｋ块在时间维上的相关性㊂而在时间维上利用ＳＶＤ与ＭＭＳＥ求Ｗ矩阵是独立的，无法直接用于算法简化，为此，本文利用时间相关性，并运用插值算法直接估计Ｗ矩阵㊂具体的，由于同ｊ不同ｋ的块在时间维上的相关性，在经过函数Ｗ＾＝ｆ（Ｈ）后，具有相关性的输入Ｈ，与得到的Ｗ之间将保持相关性㊂因此，可以利用同ｊ不同ｋ的Ｗ的相关性，通过插值算法获取某些ｋ值上的Ｗ矩阵㊂这将直接减少ＳＶＤ与ＭＭＳＥ求逆过程的计算数量㊂为了便于理解，将Ｌ维与Ｊ维（用户数）进行合并，因此Ｗ矩阵可以改写为三维矩阵㊂将矩阵按照Ｋ维展开，不同ｋ下标的矩阵可以由二维平面示意，其插值过程如图６所示㊂图６㊀Ｗ矩阵的插值示意图Ｆｉｇ．６㊀ＳｃｈｅｍａｔｉｃｄｉａｇｒａｍｏｆｉｎｔｅｒｐｏｌａｔｉｏｎｏｆＷｍａｔｒｉｘ３．１．２矩阵数据Ｗ的插值算法对于矩阵数据Ｗ的插值算法，采用ｌｉｎｅａｒ插值法㊁ｓｐｌｉｎｅ插值法与Ｐｃｈｉｐ插值法进行建模插值［２－４］㊂此处，引入 最大插值比 作为评估方法来评价插值性能，参数寻优的过程可以表示为：㊀㊀Ｒａｔｅ＝ｍａｘ１Ｒ{}ｓ．ｔ．ρｍｉｎ（Ｗ）ｍｉｎｌɪ｛１，２， ，Ｌ｝ｊɪ｛１，２， ，Ｊ｝ｋɪ｛１，２， ，Ｋ｝ρｌ，ｊ，ｋ（Ｗ）＞０．９９，（２１）式中，Ｒ表示每隔Ｒ个点进行一次插值运算㊂因此，式（２１）表示插值结果满足最小建模精度的约束下，使得插值数量越多的寻优目标㊂因此，为满足题设对于拟合后Ｗ矩阵对最小建模精度的要求，需选择合适的插值方法并计算最大插值比，为此进一步提出了最大插值比搜索方法，以评估在不同信道条件数据集下可用的插值参数，最大插值比搜索算法的详细过程在算法１中给出，其基本思想为通过给定的初始插值比Ｒａｔｅ，和选定的插值类型，不断迭代和逼近给定插值类型下的满足要求的最大插值比㊂当取得的插值比越大时，意味着Ｗ矩阵的更多部分可以通过在Ｋ维度上的相关性插值得到，不需要通过ＳＶＤ与ＭＭＳＥ求逆过程，这将大大减少计算的复杂度㊂此外，Ｄａｔａ集最大插值比的计算过程可以理解为适应信道的训练过程㊂在后续过程中，在外部信道条件未剧烈改变的情况下，不需要再次执行，因此最大插值比搜索可以在线下执行，不会影响算法的复杂度㊂算法１所提出的最大插值比搜索算法输入：训练数据集Ｄａｔａ，包含通过ＭＭＳＥ计算的标准Ｗ矩阵；初始化：初始插值比Ｒａｔｅ＝１／Ｒ，插值间隔Ｒ的初始值可以设定为Ｒ＝２㊂Ｒ的中止值可以设定为Ｒｍａｘ＝２００选定的插值类型：Ｌｉｎｅａｒ，Ｓｐｌｉｎｅ，Ｐｃｈｉｐ㊂执行过程：ＦｏｒＲ＜Ｒｍａｘ＝２００㊀ＦｏｒｉＮ＝１：Ｎ㊀㊀ＦｏｒｉＬＪ＝１：ＬＪ１．根据选定的插值比Ｒａｔｅ，确定插值点所在的横坐标序列ｘ，其中ｘ＝［１，Ｒ，２Ｒ， ，ｐＲ］Ｔ，ｐＲɤＫ＝３８４２．合并Ｊ个用户的Ｗ矩阵，使其降维为ＮˑＬＪˑＫ３．将Ｗ矩阵的第三维度中与ｘ重合的部分置零，新建为Ｗ＾，在Ｍａｔｌａｂ中可采用ｓｅｔｄｉｆｆ函数㊂置零部分准备后续进行插值填充４．执行插值操作Ｗ＾（ｉＮ，ｉＬＪ，１：Ｋ）＝ｉｎｔｅｒｐ１（ｚ，ｙ，１：Ｋ，插值类型）其中ｉｎｔｅｒｐ１表示插值函数，具体使用方法可参考Ｍａｔｌａｂ中对应函数５．评估插值结果ρｌ，ｊ，ｋ（Ｗ）＝Ｗ＾Ｈｌ，ｊ，ｋＷｌ，ｊ，ｋ ２Ｗ＾ｌ，ｊ，ｋ ２ Ｗｌ，ｊ，ｋ ２，ｌ＝１，２， ，Ｌ若ρｍｉｎ（Ｗ） ｍｉｎｌɪ｛１，２， ，Ｌ｝ｊɪ｛１，２， ，Ｊ｝ｋɪ｛１，２， ，Ｋ｝ρｌ，ｊ，ｋ（Ｗ）＞０．９９，则跳出循环（ｂｒｅａｋ）㊀㊀Ｅｎｄｉｆ㊀ＥｎｄｉｆＥｎｄｉｆ输出：插值类型，最大插值比Ｒａｔｅ表２给出了３种常见插值方法在６个Ｄａｔａ集的最大插值比㊂由于Ｄａｔａ１ Ｄａｔａ６来自不同的信道条件，因此同一插值方法在插值过程中计算出的最大插值比有较大区别㊂例如Ｄａｔａ３数据集的最大插值比显著较小，这可以理解为信道的时变性强或受到干扰噪声较大，插值算法在此时难以满足要求，需要降低最大插值比㊂进一步，横向对比３种插值方法，可见在多数的数据集下，Ｌｉｎｅａｒ插值的最大插值比最小，性能最差，这是由于简单的Ｌｉｎｅａｒ插值精度较低㊂Ｓｐｌｉｎｅ插值与Ｐｃｈｉｐ插值的最大插值比性能相近，Ｓｐｌｉｎｅ插值在Ｄａｔａ１与Ｄａｔａ３上表现比Ｐｃｈｉｐ插值较好㊂从原理上分析，可以理解为当基础函数振荡时，Ｓｐｌｉｎｅ比Ｐｃｈｉｐ更好地捕获点之间的移动，后者会在局部极值附近急剧扁平化，这在该场景下带来了更好的插值性能［５－６］㊂表２㊀３种插值方法在６个Ｄａｔａ集的最大插值比Ｔａｂ．２㊀Ｍａｘｉｍｕｍｉｎｔｅｒｐｏｌａｔｉｏｎｒａｔｉｏｏｆｔｈｅｔｈｒｅｅｉｎｔｅｒｐｏｌａｔｉｏｎｍｅｔｈｏｄｓｉｎ６ｄａｔａｓｅｔｓ最大插值比Ｌｉｎｅａｒ插值Ｓｐｌｉｎｅ插值Ｐｃｈｉｐ插值Ｄａｔａ１１／２０１／１５１／２０Ｄａｔａ２１／３０１／３０１／３０Ｄａｔａ３１／１３０１／１１１１／１１８Ｄａｔａ４１／３１／２１／２Ｄａｔａ５１／９１／９１／９Ｄａｔａ６１／１７１／１０１／１０复数矩阵运算可拆解为基本的复数运算，而基本的复数运算又可进一步拆解为基本的实数运算，实数基本运算的复杂度按照表１计算㊂对于Ｌｉｎｅａｒ插值法，根据上述描述，可得需要加减法６次㊁乘法２次㊁倒数１次，且对于复数，幅度和相位要分别插值，总复杂度是实数插值的两倍㊂然而，特殊的是这里插值点的横坐标是均匀的，因此计算的复杂度可以大大简化，仅需要３次加法㊁２次乘法和１次倒数，因此总复杂度为６８㊂每个插值时刻ｋ，共需要ＮＪＬ次插值，因此需要复杂度６８ＮＪＬ㊂本文的参数取值为Ｎ＝６４，Ｊ＝４，Ｌ＝２，因此计算复杂度为１７４０８㊂对于Ｓｐｌｉｎｅ插值法，其计算步骤为：求三次函数的系数，然后将插值点横坐标代入三次函数，计算又需要６次加法㊁６次乘法㊂且对于复数，幅度和相位要分别插值，总复杂度是实数插值的２倍，因此总复杂度为１２４，每个插值时刻ｋ，计算复杂度１２４ＮＪＬ㊂取本文参数，计算复杂度为３１７４４㊂对于Ｐｃｈｉｐ插值，由于其性能不如Ｓｐｌｉｎｅ且计算复杂度相似，因此不在此处进行考虑㊂３．２㊀降低矩阵求逆的计算复杂度对于求解逆矩阵Ｖｋ（ＶＨｋＶｋ＋σ２Ｉ）－１过程中的计算复杂度，当使用高斯消元法时［７］，求解维度ＬＪˑＬＪ的矩阵的逆矩阵的复杂度近似为Ｏ（（ＬＪ）３）；当矩阵求逆过程中使用的矩阵乘法使用文献［８］中的Ｓｔｒａｓｓｅｎᶄｓ方法时，其提出的矩阵相乘公式将常规的矩阵相乘的运算量减少很多，可以将上述复杂度降低到Ｏ（（ＬＪ）２．８０７）㊂为此，进一步研究矩阵求逆降低复杂度算法，明显看出Ｖｋ（ＶＨｋＶｋ＋σ２Ｉ）－１为Ｈｅｒｍｉｔｅ正定阵［９］，采用了一种改进的Ｓｔｒａｓｓｅｎ矩阵求逆算法［１０］，该算法结合Ｓｔｒａｓｓｅｎ矩阵求逆的高效性以及Ｈｅｒｍｉｔｅ正定阵的共轭对称性特点，使得算法运算量小，结构也简化许多㊂首先，对于矩阵分块直接求逆，假设一个Ｎ阶（这里的Ｎ被重新定义）的方阵Ａ，分块如下：Ａ＝［ａ１１］ｎˑｎ［ａ１２］ｎˑｍ［ａ２１］ｍˑｎ［ａ２２］ｍˑｎæèçöø÷ＮˑＮ㊂（２２）设Ａ的逆矩阵分块如下：Ａ＝［ｃ１１］ｎˑｎ［ｃ１２］ｎˑｍ［ｃ２１］ｍˑｎ［ｃ２２］ｍˑｎæèçöø÷ＮˑＮ，（２３）则根据矩阵分块求逆的原理有：ｃ１１＝（ａ１１－ａ１２ˑａ１２１ˑａ２１）－１ｃ１２＝－ｃ１１ˑａ１２ˑａ－１２２ｃ２１＝－ａ－１２２ˑａ２１ˑｃ１１ìîíïïï，（２４）式中，对于ＮˑＮ阶矩阵，需要的矩阵相乘的运算量级为Ｎ３㊂对于ｎ＞１，ｍ＞１，直接分块求逆算法在具体实现中需要利用递归实现，具体的运算量按照复数乘加次数统计，对于上述的直接求逆算法，以乘加次数统计理论运算量，式（２４）各部分运算量：ｃ１１的运算量为４ｍ２ｎ＋４ｍｎ２－ｍｎ，ｃ１２和ｃ２１的运算量相等，均为４ｍ２ｎ＋４ｍｎ２－６ｍｎ，同理，ｃ２２的运算量为４ｍ２ｎ＋４ｍｎ２－ｍ２，为此，矩阵Ａ利用一次分块求逆的总的运算量为：㊀㊀㊀Ｔ（１）（Ｎ）＝Ｔ（ｍ）＋Ｔ（ｎ）＋１６ｍ２ｎ＋１２ｍｎ２－１３ｍｎ－ｍ２＋４ｍ３，（２５）式中，Ｔ（１）（Ｎ）表示利用一次矩阵分块求逆算法计算矩阵求逆的总计算量，Ｔ（ｍ）和Ｔ（ｎ）分别表示对ｍ和ｎ阶复矩阵求逆所需的运算量㊂由式（２４）可知，经过一次分块求逆之后的运算量依然很高，即Ｔ（１）（Ｎ） Ｏ［（ｍａｘ（ｍ，ｎ））３］，同样可知，式（２４）中的（ａ１１－ａ１２ˑａ－１２２ˑａ２１）－１和ａ－１２２可以继续作为需要求逆的复矩阵，利用式（２２） （２４）可以继续分块求逆，所需运算量即式（２５）中Ｔ（ｍ）和Ｔ（ｎ）部分㊂采用改进的Ｓｔｒａｓｓｅｎ矩阵求逆算法，结合式（２２） （２３），Ｓｔｒａｓｓｅｎ算法应用到求逆运算有如下公式：Ｒ１＝ａ－１１１Ｒ２＝ａ２１ˑＲ１Ｒ３＝Ｒ１ˑａ１２Ｒ４＝ａ２１ˑＲ３Ｒ５＝Ｒ４－ａ２２Ｒ６＝Ｒ－１５ｃ１２＝Ｒ３ˑＲ６ｃ２１＝Ｒ６ˑＲ２Ｒ７＝Ｒ３ˑｃ２１ｃ１１＝Ｒ１－Ｒ７ｃ２２＝－Ｒ６ìîíïïïïïïïïïïïïïïïï㊀㊀㊂（２６）式中，对于ＮˑＮ阶矩阵，需要的矩阵相乘的运算量级为Ｎｌｏｇ６２＝Ｎ２．５８５，相比于矩阵分块直接求逆，运算量随着矩阵维数将有显著降低㊂按照前面相同的运算量统计方法，根据式（２５），矩阵利用一次分块求逆的总运算量为：㊀㊀㊀㊀Ｔ（１）＝Ｔ（ｍ）＋Ｔ（Ｎ）＋１３ｍ２ｎ＋１１ｍｎ２－４ｍｎ＋３ｍ２㊂（２７）对于ｎ＞１，ｍ＞１，Ｓｔｒａｓｓｅｎ矩阵求逆算法也是利用递归实现的，但因为Ｓｔｒａｓｓｅｎ算法减少了矩阵复乘次数，所以相比直接分块的常规算法运算量有明显降低㊂又由于ａ１１，ａ２２，ａ－１２２均为Ｈｅｒｍｉｔｅ矩阵，且ａＨ１２＝ａ２１，代入到式（２４）中得Ｒ３＝ＲＨ２，ｃ２１＝ｃＨ１２，根据Ｈｅｒ⁃ｍｉｔｅ矩阵的共轭对称性，式（２６）可进一步改写为：Ｒ１＝ａ－１１１Ｒ２＝ａ２１ˑＲ１Ｒ３＝ＲＨ２Ｒ４＝ａ２１ˑＲ３Ｒ５＝Ｒ４－ａ２２Ｒ６＝Ｒ－１５ｃ１２＝Ｒ３ˑＲ６ｃ２１＝ｃＨ１２Ｒ７＝ＲＨ２ˑｃ２１ｃ１１＝Ｒ１－Ｒ７ｃ２２＝－Ｒ６ìîíïïïïïïïïïïïïïïïï，（２８）式中，对于ＮˑＮ阶矩阵，需要的矩阵相乘的运算量级为Ｎ２㊂为了便于比较新求逆算法的性能改善，按照前面相同的运算量统计方法，根据式（２８），矩阵Ａ利用一次分块求逆的总的运算量为：Ｔ（１）（Ｎ）＝Ｔ（ｍ）＋Ｔ（Ｎ）＋８ｍ２ｎ＋８ｍｎ２＋ｍｎ，（２９）式中，Ｔ（１）（Ｎ）表示利用改进算法计算一次矩阵求逆的运算量，Ｔ（ｍ）和Ｔ（ｎ）部分表示对ｍ和ｎ阶复矩阵求逆所需的运算量㊂由式（２８）可知，经过一次分块求逆之后的运算量Ｔ（１）（Ｎ） Ｏ［（ｍａｘ（ｍ，ｎ））３］㊂同样可以得出，式（２８）中的Ｒ－１５可以继续作为需要求逆的复矩阵，利用式（２８）可以继续分块求逆，所需运算量即式（２９）中Ｔ（ｍ）和Ｔ（ｎ）部分㊂和矩阵直接分块求逆算法相比，新的求逆算法虽然增加了一些加减运算，但复乘次数降低㊂对于维数较高的矩阵，其中有大量的复矩阵运算，复乘消耗的运算量将远大于加减法，而这个运算量随着矩阵维数增加将有显著增大，因此新算法对于复乘次数的降低将显著改善求逆运算的实时性能㊂和常规Ｓｔｒａｓｓｅｎ矩阵求逆算法相比，改进的算法由于利用了求逆矩阵的特点，即对Ｈｅｒｍｉｔｅ矩阵进行求逆运算，所以在运算量和算法复杂度上都有明显的降低㊂常规算法计算一次矩阵求逆的计算复杂度如表３所示，改进算法计算一次矩阵求逆的计算复杂度如表４所示㊂表３㊀常规算法计算一次矩阵求逆的计算复杂度Ｔａｂ．３㊀Ｃｏｍｐｕｔａｔｉｏｎａｌｃｏｍｐｌｅｘｉｔｙｏｆｔｈｅｃｏｎｖｅｎｔｉｏｎａｌａｌｇｏｒｉｔｈｍｔｏｃａｌｃｕｌａｔｅｔｈｅｉｎｖｅｒｓｉｏｎｏｆａｍａｔｒｉｘ单项乘法次数加法次数其他Ｒ１Ｔ（ｎ）求逆Ｒ２４ｍｎ２４ｍｎ２－２ｍｎＲ３４ｍｎ２４ｍｎ２－２ｍｎＲ４４ｍ２ｎ４ｍ２ｎ－２ｍ２Ｒ５２ｍ２Ｒ６Ｔ（ｍ）求逆Ｒ７４ｍｎ２４ｍｎ２－２ｎ２ｃ１１２ｎ２ｃ２１４ｍ２ｎ４ｍ２ｎ－２ｍｎｃ１２４ｍ２ｎ４ｍ２ｎ－２ｍｎｃ２２ｍ２复杂度合计９ｍｎ（４ｍ＋４ｎ－１）＋３ｍ２＋Ｔ（ｎ）＋Ｔ（ｍ）表４㊀改进算法计算一次矩阵求逆的计算复杂度Ｔａｂ．４㊀Ｉｍｐｒｏｖｅｄａｌｇｏｒｉｔｈｍｔｏｃａｌｃｕｌａｔｅｔｈｅｃｏｍｐｕｔａｔｉｏｎａｌｃｏｍｐｌｅｘｉｔｙｏｆａｍａｔｒｉｘｉｎｖｅｒｓｉｏｎ单项乘法次数加法次数其他Ｒ１Ｔ（ｎ）求逆Ｒ２４ｍｎ２４ｍｎ２－２ｍｎＲ３Ｒ４４ｍ２ｎ４ｍ２ｎ－２ｍ２Ｒ５２ｍ２Ｒ６Ｔ（ｍ）求逆Ｒ７４ｍｎ２４ｍｎ２－２ｎ２ｃ１１２ｎ２ｃ２１ｃ１２４ｍ２ｎ４ｍ２ｎ－２ｍｎｃ２２ｍ２复杂度合计８ｍｎ４ｍ＋４ｎ－１２()＋３ｍ２＋Ｔ（ｎ）＋Ｔ（ｍ）在本文中，当矩阵维度为８（Ｊ＝４）时，改进算法总共复杂度为４７７６，常规算法总共５２５５，复杂度降低１０．０３％㊂进一步，图７给出了不同用户数量Ｊ时，改进的Ｓｔｒａｓｓｅｎ算法与常规Ｓｔｒａｓｓｅｎ算法的复杂度比较㊂可以看出，随着用户数量增加，矩阵求逆时的维度增加，采用改进的Ｓｔｒａｓｓｅｎ算法对复杂度的降低更加明显㊂图７㊀在不同用户个数Ｊ下，改进的Ｓｔｒａｓｓｅｎ算法与常规Ｓｔｒａｓｓｅｎ算法的复杂度比较Ｆｉｇ．７㊀ＵｎｄｅｒｄｉｆｆｅｒｅｎｔｎｕｍｂｅｒｏｆｕｓｅｒｓＪ，ｔｈｅｃｏｍｐｌｅｘｉｔｙｃｏｍｐａｒｉｓｏｎｏｆｔｈｅｉｍｐｒｏｖｅｄＳｔｒａｓｓｅｎａｌｇｏｒｉｔｈｍａｎｄｔｈｅｃｏｎｖｅｎｔｉｏｎａｌＳｔｒａｓｓｅｎａｌｇｏｒｉｔｈｍ３．３㊀所提算法对最小建模精度的影响利用改进的Ｓｔｒａｓｓｅｎ矩阵求逆算法求得的Ｗ＾Ｈｌ，ｊ，ｋ与原来的矩阵求逆算法求得的Ｗｌ，ｊ，ｋ进行建模估计精度计算，对于所给的数据集Ｄａｔａ１ Ｄａｔａ６，仿真不同数据集的最小建模精度，发现各数据集最小建模精度均为１，如图８所示㊂可见所提改进Ｓｔｒａｓｓｅｎ矩阵求逆算法不会对建模精度带来影响，这是因为该算法采用了分块迭代方法在变换的过程中不会带来计算误差㊂图８㊀采用改进的矩阵求逆算法后对各数据集最小建模精度的影响Ｆｉｇ．８㊀Ｉｎｆｌｕｅｎｃｅｏｆｔｈｅｉｍｐｒｏｖｅｄｍａｔｒｉｘｉｎｖｅｒｓｉｏｎａｌｇｏｒｉｔｈｍｏｎｔｈｅｍｉｎｉｍｕｍｍｏｄｅｌｉｎｇａｃｃｕｒａｃｙｏｆｅａｃｈｄａｔａｓｅｔ３．４㊀综合复杂度分析本节分别对插值算法和改进的矩阵求逆的综合复杂度进行分析㊂利用相关矩阵组的关联度降低计算复杂度，即通过Ｓｐｌｉｎｅ插值操作降低ＳＶＤ的总复杂度㊂其中，所需的乘法次数（５ＭＮ２－ＭＮ）㊁加法次数（３ＭＮ２－ＭＮ）㊁除法次数（０．５Ｎ（Ｎ－１）＋２ＭＮ）以及平方根次数（ＭＮ２）㊂最终的复杂度合计为Ｎｉｔｅ（４３ＭＮ２＋７５２Ｎ（Ｎ－１）＋１４８ＭＮ）㊂另外，通过改进的矩阵求逆，即采用改进的Ｓｔｒａｓｓｅｎᶄｓ矩阵求逆进一步降低矩阵求逆Ｖｋ（ＶＨｋＶｋ＋σ２Ｉ）－１的复杂度，得出了改进后矩阵求逆算法后所需的乘法次数（４ＭＬ２＋４Ｎ２ＬＪ）㊁加法次数（４ＭＬ２＋４Ｎ２ＬＪ－２ＭＬ－２Ｎ２）以及求逆复杂度（４７７６）㊂最终的复杂度合计为２Ｎ２（８ＬＪ－１）＋２ＭＬ（８Ｌ－１）＋４７７６㊂对于不采用插值算法的情况，每个插值时刻ｋ，由Ｈ矩阵到Ｗ矩阵，需要进行Ｊ（Ｊ＝４）次ＳＶＤ和１次ＭＭＳＥ均衡（主要复杂度在于求逆）的计算㊂其中４次ＳＶＤ需要３５７４４００次计算㊂ＭＭＳＥ均衡需要５２１１１２次计算，因此共需要计算复杂度Ｃ１＝４０９５５１２㊂采用Ｓｐｌｉｎｅ插值时，每个插值时刻的复杂度为３１７４４㊂可以计算出每次插值的复杂度收益为Ｃ２＝３５７４４００＋５２１１１２－３１７４４＝４０６３７６８㊂因此采用插值的最终复杂度收益为：ΔＣ＝（Ｃ１－Ｃ２）ˑＫˑＲａｔｅ㊂（３０）假设Ｒａｔｅ＝１／３时，ΔＣ＝５２４２２５５３６㊂可见，插值对于计算复杂度的降低比较明显㊂同样，计算复杂度可以降低为：ΔＲＣ＝Ｃ１（１－Ｒａｔｅ）＋Ｃ２ˑＲａｔｅＣ１㊂（３１）当Ｒａｔｅ＝１／３时，计算复杂度降低为原来的６６．９３％㊂当Ｒａｔｅ＝１／１０时，计算复杂度降低为原来的９０．０８％㊂综上，当采用改进的Ｓｔｒａｓｓｅｎ矩阵求逆算法时，复杂度降低了１０．０３％㊂进一步采用插值算法后，计算复杂度能够在上述的基础上再降低９．９２％（插值比为１／１０）㊁３３．０７％（插值比为１／３）㊂４㊀结论本论文主要解决ＭＩＭＯ场景下的相关矩阵组的低复杂度计算问题，首先利用Ｈ矩阵在时间相关性推导了Ｗ矩阵的相关性，通过对已有Ｗ矩阵的相关性直接插值出部分缺失的Ｗ㊂这使得在接收Ｈ矩阵时，在求取部分Ｗ矩阵后通过相关性重建完整的Ｗ矩阵；也避免了一部分Ｈ矩阵的存储以及这部分Ｈ矩阵计算ＳＶＤ与求逆获得Ｗ矩阵的过程㊂相比ＳＶＤ与求逆的复杂度，插值的复杂度要低得多㊂此外，采用了一种改进的Ｓｔｒａｓｓｅｎ矩阵求逆算法来降低求逆过程的复杂度㊂该算法结合了Ｓｔｒａｓｓｅｎ矩阵求逆的高效性以及Ｈｅｒｍｉｔｅ正定阵的共轭对称性特点，结构更简化㊂参考文献［１］㊀ＲＵＳＥＫＦ，ＰＥＲＳＳＯＮＤ，ＬＡＵＢＫ，ｅｔａｌ．ＳｃａｌｉｎｇＵｐＭＩＭＯ：ＯｐｐｏｒｔｕｎｉｔｉｅｓａｎｄＣｈａｌｌｅｎｇｅｓｗｉｔｈＶｅｒｙＬａｒｇｅＡｒｒａｙｓ［Ｊ］．ＩＥＥＥＳｉｇｎａｌＰｒｏｃｅｓｓｉｎｇＭａｇａｚｉｎｅ，２０１３（３０）１：４０－６０．［２］㊀蔡锁章，杨明，雷英杰．数值计算方法［Ｍ］．２版．北京：国防工业出版社，２０１６．［３］㊀许小勇，钟太勇．三次样条插值函数的构造与Ｍａｔｌａｂ实现［Ｊ］．兵工自动化，２００６（１１）：７６－７８．［４］㊀陈帅，岳迎春，徐巍，等．小波时间序列对非平稳信号中突变点的辨识与处理［Ｊ］．测绘科学，２０１３，３８（５）：１１－１２．［５］㊀ＤＥＢＯＯＲＣ．ＡＰｒａｃｔｉｃａｌＧｕｉｄｅｔｏＳｐｌｉｎｅｓ［Ｊ］．ＡｐｐｌｉｅｄＭａｔｈｅｍａｔｉｃａｌＳｃｉｅｎｃｅｓＮｅｗＹｏｒｋＳｐｒｉｎｇｅｒ，１９７８，２７（１４９）：１５７－１５７．［６］㊀ＦＲＩＴＳＣＨＦＮ，ＣＡＲＬＳＯＮＲＥ．ＭｏｎｏｔｏｎｅＰｉｅｃｅｗｉｓｅＣｕｂｉｃＩｎｔｅｒｐｏｌａｔｉｏｎ［Ｊ］．ＳｉａｍＪｏｕｒｎａｌｏｎＮｕｍｅｒｉｃａｌＡｎａｌｙｓｉｓ，１９８０，１７（２）：２３８－２４６．［７］㊀颜志升，郑昱．基于高斯消元的自适应信号处理的实现方法：ＣＮ１１１４２７０１４Ａ［Ｐ］．２０２０－０７－１７．［８］㊀ＳＴＲＡＳＳＥＮＶ．ＧａｕｓｓｉａｎＥｌｉｍｉｎａｔｉｏｎｉｓＮｏｔＯｐｔｉｍａｌ［Ｊ］．ＮｕｍｅｒｉｓｃｈｅＭａｔｈｅｍａｔｉｋ，１９６９，１３（４）：３５４－３５６．［９］㊀杨忠鹏，林志兴．关于Ｈｅｒｍｉｔｉａｎ矩阵的特征的注记［Ｊ］．大学数学，２００３，１９（５）：５２－５３．［１０］李瑞，李晓明，董晔．ＳＴＡＰ中的协方差矩阵求逆快速算法研究［Ｊ］．计算机仿真，２０１１，２８（２）：２５－２８．作者简介：㊀㊀施育鑫㊀国防科技大学第六十三研究所博士研究生㊂主要研究方向：通信抗干扰㊁ＯＦＤＭ㊂㊀㊀鲁信金㊀国防科技大学电子科学学院博士研究生㊂主要研究方向：信息论㊁索引调制㊁ｐｏｌａｒ码㊁物理层安全㊁无线通信技术等㊂在各类期刊和会议论文集上发表论文多篇㊂㊀㊀孙艺夫㊀国防科技大学电子科学学院博士研究生㊂主要研究方向：可重构信息超表面㊁通信抗干扰㊁物理层安全㊂㊀㊀雷㊀菁㊀国防科技大学电子科学学院教授，博士生导师㊂主要研究方向：信息论㊁ＬＤＰＣ码㊁物理层安全㊁ｐｏｌａｒ码㊁无线通信等㊂㊀㊀李玉生㊀国防科技大学第六十三研究所正高级工程师，硕士生导师㊂主要研究方向：通信抗干扰㊂。

基于稀疏表示的频控阵 MIMO 雷达多目标定位陈慧;邵怀宗;潘晔;王文钦【摘要】针对频控阵多输入多输出（MIMO）雷达，提出了一种基于压缩感知稀疏表示思想的目标定位算法。

首先回顾了 MIMO 雷达和频控阵的特点，进而研究了频控阵 MIMO 雷达的性质，它不但可以具有MIMO 雷达的优点，而且能够感知目标的距离维信息，同时针对频控阵 MIMO 雷达接收数据模型进行数学建模，并把目标定位问题表示成稀疏表示框架下的代价函数。

最后利用凸优化工具对代价函数进行优化求解，由所得稀疏权向量中的非零元素索引映射出目标的方位和距离信息。

与现有的经典 MUSIC 算法相比，具有更好的定位性能，计算机仿真结果证明了所提算法的有效性。

%For the frequency diverse array(FDA)multiple-input and multiple-output(MIMO)radar,a target localization algorithm in sparse signal representation perspective is presented.Firstly,the characters of the MIMO radar and FDA are reviewed,then the properties of the FDA MIMO radar are studied,that is,it not only owns the merits of the MIMO radar,but also can sense the range information of the targets.Its re-ceiving mathematical measurement is also modeled,and the target localization problem is described as a cost function under the sparse representation framework.Finally,the angle and the range of the targets are esti-mated by mapping the non-zero element indexes of the sparse vector which is obtained by solving the cost function using existing convex pared with the existing classic MUSIC algorithm,the proposed algo-rithm can achieve better localization performance.The computer simulation results demonstrate the effective-ness of the algorithm.【期刊名称】《雷达科学与技术》【年(卷),期】2015(000)003【总页数】6页(P259-264)【关键词】频控阵 MIMO 雷达;压缩感知;稀疏表示;参数估计【作者】陈慧;邵怀宗;潘晔;王文钦【作者单位】电子科技大学通信与信息工程学院，四川成都 611731;电子科技大学通信与信息工程学院，四川成都 611731;电子科技大学通信与信息工程学院，四川成都 611731;电子科技大学通信与信息工程学院，四川成都 611731【正文语种】中文【中图分类】TN9580 引言多输入多输出(Multiple-Input and Multiple-Output,MIMO)雷达近年来受到了广泛关注[1-3]。

多输入多输出线性定常系统稳定裕度的分析与改进李信栋;苟兴宇【摘要】针对多输入多输出(MIMO)控制系统的稳定裕度求解问题,首先分析了现有的回差阵奇异值法这一计算方法,并得到其解决单输入单输出(SISO)系统的稳定裕度结论,在此基础上,提出两种基于系统回差阵的稳定裕度改进方法;一种是在有限条件下利用矩阵的特征值代替奇异值来建立与稳定裕度关系的策略,另一种是利用系统逆回差阵的行列式,通过求其奇异值来计算系统稳定裕度;最后结合工程实例,通过数值仿真验证两种稳定裕度计算方法相比原方法都有不同程度的改进,而且三种方法可以结合起来进行分析,最大化的减小系统稳定裕度结果的保守性.【期刊名称】《控制理论与应用》【年(卷),期】2014(031)001【总页数】7页(P105-111)【关键词】回差阵奇异值法;特征值;逆回差阵;稳定裕度【作者】李信栋;苟兴宇【作者单位】北京控制工程研究所,北京100190;空间智能控制技术重点实验室,北京100190【正文语种】中文【中图分类】TP273稳定裕度是判断系统相对稳定性的一个重要依据.在经典控制理论中,单输入单输出(single-input singleoutput,SISO)线性定常系统的稳定裕度常用频域下的相角裕度和幅值裕度来度量,并且有明确的计算方法;在工程中,控制器参数的调节与优化是一个很重要的问题,文献[1]就根据改进的D--分割法及与最大灵敏度、超调和调节时间有关的目标函数,给出一种基于幅值裕度和相位裕度的比例--积分--微分(proportionintegral-derivative,PID)参数最优整定方法.SISO系统稳定裕度已能有效地应用于工程实际.然而,对于多输入多输出(multi-input multi-output,MIMO)线性系统由于其各回路之间存在耦合,以及频域理论的应用限制,目前没有公认明确的稳定裕度结论.最近几十年众多学者另辟蹊径,基于现代控制技术对MIMO系统稳定裕度的研究做了大量工作.利用频域下的Nyquist理论,Latchman和Crisalle等人提出临界方向理论(critical direction theory)[2–3],通过定义临界方向和临界扰动半径等概念来计算系统Nyquist鲁棒稳定裕度,其中文献[3]将这一理论推广至MIMO线性系统.Doyle通过定义结构奇异值[4]µ(ω)来分析结构不确定性MIMO线性系统稳定裕度;文献[5]依据µ分析法处理多回路不确定性问题所具有的独特优势,建立了µ(ω)值与经典增益、相位裕度之间的转换关系.纪多红等人[6]提出了基于µ分析的多变量系统稳定裕度评估方法,评估结果的精确性取决于所估计µ值的精确程度,若估计的µ值不准确,所得稳定裕度必然具有较大保守性.与结构奇异值µ(ω)分析法类似,de Gaston R R E和Safonov G等人[7–8]提出一种实多回路稳定裕度计算方法,基于映射理论通过迭代算法得到MIMO系统稳定裕度结果km,针对实参数扰动可写成对角矩阵的情形,文献[7]将稳定裕度的计算问题转化为计算km的上下界问题;文献[8]针对一类实参不确定性系统,解决了不确定性元素之间相关联时系统的稳定裕度问题,虽然实多回路稳定裕度理论可以得到比较准确的稳定裕度结果,但是同时也可看出这种方法的复杂程度,相对于其它MIMO系统稳定裕度的求取方法,其要求的计算量非常大.相比较而言,文献[9–10]提出的基于回差阵奇异值计算稳定裕度的策略所要求的计算量较小,结合奇异值理论,通过分析系统回差阵即可得到系统的稳定裕度,其中文献[10]将所提方法应用到某飞行器的两输入两输出横侧向姿态控制系统.文献[11]在回差阵奇异值法基础上进行了改进,同样应用到某型飞机的侧向运动系统进行了稳定裕度的分析.鉴于回差阵奇异值法在计算MIMO系统稳定裕度方面的优越性,本文将对这一方法进行详细分析,在此方法的基础上作出改进,以减小系统稳定裕度计算的保守性.首先将这一理论退化为解决SISO线性定常系统的情形,得到相应的稳定裕度结论;其次尝试用矩阵的特征值代替奇异值来分析系统回差阵,得到基于回差阵特征值的稳定裕度改进方案;最后根据逆回差阵信息,结合矩阵奇异值理论得到MIMO系统稳定裕度计算的另一种改进方法.将3种方法所得稳定裕度结合起来分析,最大化的减小系统稳定裕度结果的保守性.奇异值分解是矩阵分解中的一个重要的技术手段,是现代数值分析的最基本和最重要的工具之一,可为控制系统的分析和设计提供有力支持.下面对回差阵奇异值法[9–10]这一重要的稳定裕度理论进行分析.系统模型如图1所示,其中G(s)∈ℂn×n为系统标称模型,则系统回差阵为I+G(s).在输入端引入不确定性指数阵形式量测阵可知系统回差阵变为I+G(s)P(s),下面分析当所有回路中的ki或φi同时变化时,系统仍能保证自身稳定性所容许其变化的最大值.已知对于标称模型的回差阵[12],有其中:ϕc(s)是闭环特征多项式,ϕo(s)是开环特征多项式,α是常数.若闭环标称系统是渐近稳定的,说明det[I+G(s)]的分子ϕc(s)的根全都具有负实部,因此令s=jω,必有det[I+G(s)]/=0,根据奇异值的性质,则意味着即回差矩阵在工作频率范围内的最小奇异值就是稳定裕度的衡量标准.同理,若引入量测阵后的系统也是渐近稳定的,则有本文研究的重点是系统由渐近稳定达到临界稳定时,量测阵P(s)中ki和φi所能变化的最大范围.根据P(s)的定义知ki不可能等于零,因此P是非奇异的,而又已知I+G非奇异,利用矩阵分离特性可得进一步可得根据奇异值的性质分析式(4),可知若成立,则可得系统稳定的充分条件为若考虑增益和相位在每个通道同时变化的情况,应用方程(1),可得确保MIMO线性定常系统稳定的充分条件将回归矩阵最小奇异值作为参数,若G)≥m,根据式(7),分别令φ=0◦和k=1,可以得到多变量系统每个通道同时确保的幅值裕度GM和相位裕度PM表达式:通过上面的分析可看出,保证系统稳定的条件式(4)是一个保守条件,因此不可避免地,采用回差阵奇异值法来求取系统的稳定裕度也会导致保守的结果,这是此方法仍需改进之处.2.1 退化结果(Degradation results)本节分析MIMO线性系统的稳定裕度结果退化为SISO系统的情形.对于SISO线性系统g(s),不确定性结构相应的由不确定性量测阵变为如下标量形式:接下来探讨系统达到临界稳定状态时p(s)中参数k和φ所能变化的最大范围.首先,笔者假定标称状态下系统是稳定的,即满足条件引入不确定性p(s),若使系统仍保持稳定,则有利用分离特性:根据p(s)的定义可知p(s)＞0,结合式(10),由上式可得分析知上述不等式等价于对第一个条件进行整理，可等价为或考虑到系统g(jω)的一般性,根据此条件不能得到明确的稳定裕度量测值p(jω)使得式(11)成立.类似地,对于第2个条件式,可整理得上式说明稳定裕度量测值p(jω)随g(jω)变化,同样不能得到满足式(11)的明确稳定裕度结果.因此,对于一般系统g(jω),可知使式(11)成立的充分条件为与上文中MIMO线性系统分析类似,这里可得系统稳定的一个充分条件:考虑回路中增益和相位的变化,将方程(9)代入上式,可得确保系统稳定的充分条件: 将min[1+g(jω)]作为参数,若min[1+g(jω)]≥m,根据上式,则可得SISO系统确保的幅值裕度GM和相位裕度PM为从上式所示结论可看出,回差阵奇异值法是在一个频率点处得到了系统增益裕度和相位裕度两者的计算表达式,这和经典控制里SISO系统稳定裕度的定义及计算结果不同,经典控制理论中需要两个频率点(即截止频率和穿越频率)来分别确定系统的相位裕度和增益裕度两个裕度指标,可认为此处差异是回差阵奇异值法计算结果的保守性在SISO线性定常系统下的具体体现.然而同样基于这一点可以进行如下大胆推测:对于MIMO线性定常系统,可否通过定义一个指标即可有效地衡量整个系统的相对稳定程度.特征值分解是分析矩阵特性的另外一个重要工具,矩阵特征值和奇异值之间有着紧密联系,其中重要一点就是特征值的幅值是最小奇异值的上界;本节在文献[10–11]的基础上,针对前面提出的回差阵奇异值法具有的保守性,下面尝试利用特征值代替奇异值来进行MIMO线性系统稳定裕度的分析,并给出详细分析过程.对于如图1所示系统,同样在输入端引入不确定性指数阵形式量测阵P(s),此时标称状态下系统是稳定的充要条件变为即利用回差矩阵在工作频率范围内的最小特征值衡量稳定裕度.若加入不确定性阵P(s)后系统仍稳定,则有利用如式(3)所示矩阵分离特性,同样假定P(s)是非奇异的,且已知I+G非奇异.因此矩阵P(s)和I+G没有为零的特征值,根据特征值与矩阵行列式的关系,可知只有下面的条件式:成立,才能够满足式(15)的不等式要求.这里需要注意的是,矩阵特征值不像奇异值那样具有类似如式(5)所示不等式的性质.因此,对于上面式(16),只有当对角阵P(s)的所有对角元素相同,即满足才有对于矩阵特征值来说,条件恒成立.在此条件下对上面式(18)进一步整理,可得因此要式(15)成立,需满足条件即上式就是使系统稳定的最终的充分条件;考虑增益和相位在每个通道同时变化的情况,根据方程(20),可得不确定系统稳定的充分条件为其中min|λ(I+G)|与稳定裕度的关系如式(8)所示.由于特征值法本身性质的不足,导致稳定裕度的计算不精确甚至是错误的.性质的不足体现在:在推导出式(20)所示系统稳定充分条件的过程中,不能使用类似于如式(5)所示不等式的奇异值性质,这种性质在特征值理论中并不是恒成立的;为推导得式(20)必须要求不确定对角阵P(s)的所有对角元素相同,即满足式(17).与文献[11]相比,本文不但给出了利用特征值代替奇异值来进行MIMO线性系统稳定裕度的分析的详细过程,并且用min|λ(I+G)|代替文献[11]中的给出了最终系统稳定充分条件,使结果更加准确.注1虽然本节在分析过程中对量测阵的形式有一定的限制条件,但是正如在开头所述,特征值的幅值是最小奇异值的上界,应用特征值幅值的最小值代替最小奇异值可明显改善系统稳定裕度结果的保守性,这一点使得本节分析结果更具有实际工程价值和意义.上文中讲到可通过计算回差阵奇异值来求取MIMO系统稳定裕度的算法,算法指出若标称系统闭环稳定,则闭环特征多项式ϕc(s)的根全具有负实部,而且根据式(2),可进一步得det[I+G(s)]/=0,这一点是采用奇异值理论分析MIMO系统稳定裕度的理论基础.根据此出发点,下面从另一个角度分析系统的稳定性及稳定裕度.首先假设标称系统是可逆的,即G−1(s)存在,系统回差阵可进行如下分解:对上式两边分别取行列式,可得结合式(2),可进一步得其中ϕz(s)是开环系统零点特征多项式.通过上式可以看出,逆回差阵I+G−1(s)同样可以与闭环特征多项式建立联系,应用奇异值理论分析det[I+G−1(s)],亦能得到系统稳定裕度,称这种方法为逆回差阵奇异值法.类似的,若标称系统闭环稳定,则有可知I+G−1(s)是可逆的,即加入不确定量测阵P(s)后系统仍保持稳定的条件为即.利用矩阵分离特性根据P(s)的定义,可知P−1(s)是非奇异的,又知I+ G−1(s)亦是非奇异的,可得式(24)成立的一个充分条件为由如式(5)所示奇异值性质,可得联系式(25),令不等式(26)右边项小于1,则不确定系统稳定的充分条件变为考虑增益和相位在每个通道同时变化的情况,结合P(s)的定义式,可得系统稳定的充分条件将回归矩阵最小奇异值作为参数,根据条件式(28),可确定系统保持稳定所容许的每个通道增益和相位的同时变化范围.具体地,分别令φ=0◦和k=1,得到MIMO线性系统允许各回路同时变化的最大幅值裕度GM和相位裕度PM分别为通过上述分析可以看出,利用系统逆回差阵同样可以得到系统的稳定裕度,将其作为回差阵奇异值法的有效补充,可减小稳定裕度计算结果的保守性.至此分析得到了基于系统回差阵及逆回差阵信息的两种稳定裕度改进算法.下面通过工程实例对文中所提出的两种MIMO线性定常系统稳定裕度改进方法—回差阵特征值法和逆回差阵奇异值法进行验证,证明两者在回差阵奇异值法的基础上都取得了不同程度的改进.例1以文献[13]中本体与天线同时机动的两刚体卫星系统为例,对所提出的新方法进行验证.针对文献中的非线性系统,在天线指向角在俯仰平面内小角度机动的假设条件下,可得其动力学方程,具体形式如下:其中:θ为卫星姿态转角,β为天线转角,Mcy为绕航天器Y轴的俯仰姿态控制力矩,βf为天线指令输入转角.若选择输入和输出分别为可得对象传递函数矩阵为设计控制器在计算稳定裕度之前,首先要判断系统的稳定性,只有稳定的系统的稳定裕度才有意义.具体地,可通过下面的定理证明系统的稳定性.定理1多变量反馈系统为渐近稳定的充分必要条件是[12]其中:n0是系统开环不稳定极点个数,enc det[I+ G(s)]表示det[I+G(s)]的Nyquist 图以顺时针方向包围原点的次数,定理指出系统稳定的充分必要条件是det[I+G(s)]的Nyquist曲线逆时针包围原点的次数等于系统的开环不稳定极点数.如果开环系统是稳定的,即n0=0,那么MIMO闭环系统稳定的充分必要条件是det[I+G(s)]的Nyquist曲线不包围原点.计算上述系统回差阵I2+G(s)K(s)及det[I2+ G(s)K(s)].明显地,系统没有右半平面开环极点,可知系统是开环稳定的.进一步作det[I2+G(s)K(s)]的Nyquist曲线如图2所示.由图2可知det[I2+G(s)K(s)]的Nyquist曲线不包围原点,说明闭环系统稳定.若设计卫星本体俯仰轴姿态期望跟踪曲线和卫星天线指向角期望跟踪曲线分别为在零初始条件下得到系统的跟踪误差曲线见图3.通过图3的仿真曲线可以看出系统是稳定的且误差趋近于0,接下来分别计算3种方法所求得的稳定裕度结果.首先采用回差阵奇异值法,根据其回差阵I2+ G(s)K(s),通过MATLAB仿真工具可求得回差阵最小奇异值为minσ(I+GK)=0.6069,根据式(8)可得相应的稳定裕度为:当所有通道的幅值不变时,允许所有通道的相位变化是35.3◦,当所有通道的相位不变时,允许所有通道的幅值变化是8.11dB.对于回差阵特征值法,同样根据回差阵I2+ G(s)K(s),然后通过仿真可求得回差阵特征值幅值的最小值为min|λ(I+GK)|=0.6109,可得相应的稳定裕度为:当所有通道的幅值不变时,允许所有通道的相位变化是35.6◦,当所有通道的相位不变时,允许所有通道的幅值变化是8.2dB.最后,采用逆回差阵奇异值法,需计算系统逆回差阵I2+[G(s)K(s)]−1,然后通过仿真工具可求得逆回差阵最小奇异值为minσ[I+(GK)−1]=0.7590,根据式(29)可得相应的稳定裕度为:当所有通道的幅值不变时,允许所有通道的相位裕度是44.6◦,当所有通道的相位不变时,允许所有通道的幅值裕度是4.9dB.为方便对3种方法所求得的稳定裕度结果进行比较,将3者的结果见表1.根据表1所示,对比稳定裕度结果可看出,两种改进后方法的相位裕度比回差阵奇异值法都有所提高,其中逆回差阵奇异值法求得的相位裕度结果最大,为44.6◦.回差阵特征值法求得的增益裕度结果最大,为8.2dB.但同时也应该看到,逆回差阵奇异值法求得的增益裕度结果较小;因此,可以将3种方法所求得的稳定裕度结果结合起来一起考虑,得到最终保守性较小的系统稳定裕度.由本文中例子中可知,系统保持稳定所容许的两个回路的增益可同时增大8.2dB,或相位同时滞后44.6◦.下面进一步对所求稳定裕度结果进行分析,通过仿真验证其精确性和有效性.根据计算结果知系统增益裕度为8.2dB,即各回路增益可同时增大2.57倍.对系统进行仿真可得跟踪误差曲线如图4所示.从图4中可以看出系统仍是稳定的,若继续增大各回路增益至2.59倍,可得系统仿真曲线如图5所示.从图5中可以看出,天线指向角误差曲线已有发散趋势,说明系统临界不稳定.通过上述仿真分析可知,利用文中改进方法所求得稳定裕度已经非常接近系统真实稳定裕度,有效地减小了计算结果的保守性.本文针对回差阵奇异值法这一计算MIMO系统稳定裕度的重要理论方法进行了详细分析,讨论这一理论退化为解决SISO线性定常系统的情形,得到相应的稳定裕度结论;然后基于回差阵信息对系统稳定裕度分析方法进行改进.首先尝试使用矩阵的特征值代替奇异值来分析系统回差阵,在某些条件限制下得到了基于回差阵特征值的稳定裕度结果;此外,利用系统逆回差阵理论,借助矩阵奇异值计算系统的稳定裕度;最后通过实例分析指出三种策略相结合使用可以有效减小所计算系统稳定裕度的保守性.本文虽然基于系统的回差阵及逆回差阵信息,针对MIMO线性定常系统稳定裕度的分析做了一定的工作,但是这些算法仍有不足;因此下一步的工作重点仍将是提出公认的更为有效的保守性更小的稳定裕度指标和算法;或者如文中所推断,设法通过定义一个新的指标来确切衡量系统的稳定裕度.李信栋(1986–),男,博士研究生,研究方向为航天器姿态控制, E-mail:*****************;【相关文献】[1]欧林林,顾诞英,张卫东.基于幅值裕度和相位裕度的PID参数最优整定方法[J].控制理论与应用,2007,24(5):837–840. (OU Linlin,GU Danying,ZHANG Weidong.Optimal tuning method of PID controller based on gain margin and phase margin[J].ControlTheory&Applications,2007,24(5):837–840.)[2]AL-SHAMALI S A,JI B W,CRISALLE O D,et al.The nyquist robust sensitivity margin foruncertain closed-loop systems[J].International Journal of Robust and Nonlinear Control,2005,15(14):619–634.[3]SON J E,MANI A S,LATCHMAN H A.Robustness analysis for MIMO aystems with unstructured uncertainties[C]//Proceedings of the 8th IEEE International Conference on Control and Automation. New York:IEEE,2010:1333–1337.[4]DOYLE J.Analysis of feedback systems with structured uncertainties[J].IEE Proceedings Part D,1982,129(6):242–250.[5]TSAO T T,LEE F C,AUGENSTEIN D.Relationship between robustnessµ-analysis and classical stability margins[C]//Proceedings of IEEE Aerospace Conference.NewYork:IEEE,1998,4:481–486.[6]纪多红,刘林,唐强.µ分析在飞控系统稳定裕度评估中的应用[J].飞行力学,2007,25(4):64–68. (JIduohong,LIUlin,TANGqiang.Applicationofµanalysisinclearance of stability margin of fl ight control system[J].Flight Dynamics, 2007,25(4):64–68.)[7]DE GASTON R R E,SAFONOV M G.Exact calculation of the multiloop stabilitymargin[J].IEEE Transactions on Automatic Control, 1988,33(2):156–171.[8]SIDERIS A,DE GASTON R R E.Multivariable stability margin calculationwithuncertaincorrelatedparameters[C]//Proceedingsofthe 25th Conference on Decision and Control.New York:IEEE,1986, 766–771.[9]LEHTOMAKI N A,SANDELL N R,ATHANS M.Robustness results in linear-quadratic gaussian based multivariable control designs[J].IEEE Transactions on Automatic Control,1981,26(1):75–92.[10]MUKHOPADHYAY V,NEWSOM J R.Application of matrix singular value properties for evaluating gain and phase margins of multiloop systems[C]//AIAA Guidance Navigation and Control Conference.California:AIAA,1982,420–428.[11]吴斌,程鹏.多变量飞控系统的稳定裕度分析[J].航空学报,1998, 19(6):657–661. (WUBin,CHENG Peng.Stability margin analysis of the multiloop fl ight control systems[J].Acta Aeronautica ET Astronautica Sinica, 1998,19(6):657–661.)[12]白方周,庞国仲.多变量频域理论与设计技术[M].北京:国防工业出版社,1988. (BAI Fangzhou,PANG Guozhong.Multivariable Frequency Domain Theory and Design Technology[M].Beijing:National Defence Industry Publishing Company,1988.)[13]曾海波.挠性多体卫星指向控制设计研究[D].北京:北京控制工程研究所,2003 (ZENG Haibo.Study on the design of fl exible multibody satellite pointingcontrol[D].Beijing:Beijing Institute of Control Engineering,2003.)。