CN104914451A - 一种块Toeplitz矩阵低复杂度求逆的空时抗干扰方法 - Google Patents

Abstract

本发明属于阵列信号处理技术领域，涉及一种块Toeplitz矩阵低复杂度求逆的空时抗干扰方法，具体步骤为：(S1)初始化，将空时权值求解问题表述为线性方程组求解问题；(S2)计算采样协方差矩阵的近似值，利用平稳性假设，保证采样协方差矩阵的估计值为块Toeplitz矩阵；(S3)以分块矩阵的形式，重新表述空时权值求解问题；(S4)定义迭代求解空时权值所需的变量和算符；(S5)将空时权值的求解分解为跌带求一系列的方程组，最终得到空时权值；本发明提出了迭代求解空时权值的低复杂度方法。与传统的SMI方法相比，在采用相同抗干扰准则的前提下，二者的性能相同，本发明方法的计算复杂度大大降低了传统方法的计算复杂度。

Description

一种块Toeplitz矩阵低复杂度求逆的空时抗干扰方法

技术领域

本发明属于阵列信号处理技术领域，具体涉及一种块Toeplitz矩阵低复杂度求逆的空时抗干扰方法，适应于抗干扰天线阵的空时滤波处理方法，特别适应于在干扰环境下工作的卫星导航接收机，也可运用在其他类型自适应阵列信号处理中。

背景技术

卫星导航信号到达地面时非常微弱，极易受到有意或无意干扰。在军事应用和航空应用等有高可靠性要求的领域，卫星导航接收机通常采用自适应天线阵技术来提高其抗干扰能力。为同时利用空域与时域自由度，抗干扰天线阵通常采用空时自适应处理(Space-time adaptive processing，简称STAP)，即采用空时自适应滤波器同时进行空域和时域滤波。在本发明中，空时滤波器所对应的阵元个数用N表示，每个阵元之后的时域抽头个数相同，用K表示。为使空时滤波器可通过自均衡的方式补偿各天线阵元及射频通道的群时延差异，以达到更好的干扰抑制和增强有用信号的效果，需要适当地增加时域抽头个数K。从理论上看，求解空时滤波器中对应各抽头的加权值(下文简称空时权值)首先需要对各抽头处数据的协方差矩阵进行估计以得到采样协方差矩阵，然后求解涉及采样协方差矩阵求逆的线性方程组。为避免直接进行矩阵求逆运算，早期的许多文献从准最优的角度提出了以多级维纳滤波为代表的降维、降秩方法。

近年来，随着数字信号处理实现技术的发展，直接快速实现矩阵求逆运算变为可能。为满足更严苛的抗干扰要求，使用直接采样矩阵求逆方法(Sample Matrix Inversion，简称SMI)的STAP得到了广泛的运用。由于采样协方差矩阵为NK×NK维，故通常的SMI方法计算量为O[N³ K³]。由此可知增加时域抽头个数，将使得STAP处理的计算量及相应的工程实现代价急剧增大，限制了其在机载、弹载等小型化武器平台上的广泛应用。针对SMI方法计算量大的问题，大量文献对减少SMI方法中空时权值求解的计算量进行了研究，现有技术主要包括：1)针对无线通信应用环境下采样协方差矩阵的特点，加速矩阵求逆的算法(参见文献“Zhu H，Chen W，SheF.Improved fast recursive algorithms for V-BLAST and G-STBC with novelefficient matrix inversion，Dresden，Germany，2009[C].Institute of Electricaland Electronics Engineers Inc，2009.”)；2)利用雷达阵列信号的STAP处理中采样协方差矩阵的Hermite对称性减少计算量(参见文献“高飞，王永良，陈辉，等.STAP中的矩阵求逆问题研究[J].雷达科学与技术，2008，6(3)：215-218.”与“Yang X，Liu Y，Long T.Pulse-order recursive methodfor inverse covariance matrix computation applied to space-time adaptiveprocessing[J].Science China Information Sciences，2013，56(4)：1-12.”)。然而，上述传统的矩阵求逆方法的计算量仍然与K³成正比，其计算量仍然会快速随时域抽头数的增加而增长。考虑到通常的STAP处理中，阵元数N较小，因此减少空时权值求解的计算量随时域抽头个数K增长的速度十分必要。

发明内容

为减少STAP处理的计算量及其工程实现代价，本发明提出了利用采样协方差矩阵可表示为块Toeplitz矩阵的特点，快速实现空时权值的求解。Toeplitz矩阵，即托普利兹矩阵，简称为T型矩阵，T型矩阵的特点是：除第一列外，其他每个元素都与左上角的元素相同，块Toeplitz矩阵是对Toeplitz矩阵概念的推广。对于一个矩阵，如果将其表示为分块矩阵的形式后，其各子块矩阵之间满足类似与Toeplitz矩阵的斜对角相等的关系，即为块Toeplitz矩阵。本发明的具体技术方案步骤如下：

(S1)初始化，将空时权值求解问题表述为如下的线性方程组求解问题：

设第a个阵元后第b个时域抽头处的复基带信号为采样值序列x_ab[j](j为整数，代表采样时刻)，空时滤波器中与之相应的加权值的复共轭为w_ab。将同一时刻各抽头处的采样值先按照阵元先后顺序排列，再按照抽头先后顺序排列，可构成如下的随机向量：

X＝[x₁₁，x₂₁，…，x_N1，x₁₂，x₂₂，…，x_N2，…，x_1K，x_2K，…，x_NK]^T    (1)

其中，N表示阵元个数，K表示抽头个数；a的取值为1，…，N之间的整数，b的取值为1，…，K之间的整数；

相应地，将所有的空时权值的复共轭可写为如下的向量形式：

W＝[w₁₁，w₁₂，…，w_1K，w₂₁，w₂₂，…，w_2K，…，w_N1，w_N2，…，w_NK]^T    (2)

根据上述定义，空时滤波器的输出可表示为：

y＝W^HX    (3)

其中，T表示矩阵转置，(·)^H表示对矩阵或向量进行共轭转置操作，下文中都采用此符号。

阵列各抽头处数据的协方差矩阵为E(XX^H)，E(·)表示求数学期望运算。工程实践中通常将数学期望用求平均值取代，由此得到作为协方差矩阵有效近似的采样协方差矩阵R_xx。设在一个权值更新周期内有L个采样点，则R_xx的第u行第v列元素表示如下：

${\left(R_{\mathrm{xx}}\right)}_{\mathrm{uv}}=\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{\left(X\right)}_u{\left(X\right)}_v^{*}---\left(4\right)$

上式中：(X)_u表示向量X中的第u个元素，(X)_v表示向量X中的第v个元素，(·)^*表示求复数的转置。

根据传统的空时抗干扰理论，按照一定的优化准则所寻找到的最优空时权值，通常可表示为如下形式的线性方程组的解。

R_xxW＝c    (5)

上式中c为由抗干扰优化准则所决定的约束向量。

(S2)将采样协方差矩阵表示为块Toeplitz矩阵的形式，简化其计算：

假设在一个权值更新周期内，各阵元上的信号是平稳的，且各阵元之间是联合平稳的，则仅取决于a₁、a₂与b₁-b₂。由此可推导出协方差矩阵可表示为块Toeplitz矩阵形式，相应地采样协方差矩阵也可以采用如下的块Toeplitz矩阵形式来计算：

$R_{\mathrm{xx}}=\left[\begin{array}{ccccc}{R}_0& R_1& R_2& ...& R_{K-1}\\ R_{-1}& R_0& R_1& ...& R_{K-2}\\ R_{-2}& R_{-1}& R_0& ...& R_{K-3}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ R_{-K+1}& R_{-K+2}& R_{-K+3}& ...& R_0\end{array}\right]---\left(6\right)$

上式中R_i，i＝-K+1，…，0，…，K-1为2K-1个N×N方阵。

易知R_xx为Hermite矩阵(Hermite矩阵，即厄米矩阵，其构成元素关于主对角线共轭对称)，因此其中的各子块之间满足如下的共轭对称性：

$R_i^H=R_{-i},i=-K+1,...,0,...,K-1---\left(7\right)$

由此可知，计算采样协方差矩阵可以简化为计算K个N×N的子阵R_m，m＝0，…，K-1。R_m的第u行第v列元素按照下式计算：

${\left(R_m\right)}_{\mathrm{uv}}=\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{x}_{u1}\left[j\right]x_{v(m+1)}^{*}\left[j\right]---\left(8\right)$

(S3)将式(5)的方程进行简化处理，使其适合于用迭代方式求解：

将空时权值的转置及约束向量的转置均表示为K个1×N的行向量构成的分块矩阵，即：

W^H＝[A₁，A₂，…，A_K]c^H＝[C₁，C₂，…，C_K]    (9)

利用R_xx的共轭对称性，将式(5)的空时权值求解问题表示为：

[C₁，C₂，…，C_K]＝[A₁，A₂，…，A_K]R_xx    (10)

(S4)定义表述迭代求解方程式(10)所需的变量及算符：

设迭代步数为p＝1，2，…，K，L_p+1为第一行子块为R₀、R_-1至R_-p，第一列子块为R₀、R₁至R_p所构成的块Toeplitz矩阵，即：

$L_{p+1}=\left[\begin{array}{ccccc}{R}_0& R_{-1}& R_{-2}& ...& R_{-p}\\ R_1& R_0& R_{-1}& ...& R_{-p+1}\\ R_2& R_1& R_0& ...& R_{-p+2}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ R_p& R_{p-1}& R_{p-2}& ...& R_0\end{array}\right]---\left(11\right)$

对于构成元素均为N×N矩阵的任意分块矩阵D，引入两个算符：一个为表示对矩阵进行以子块为单位的转置操作，即中第(u，v)子块，正好是D中的第(v，u)子块；另一个为表示交换矩阵中同一行中各子块的先后顺序，即若设： $\tilde{D}=[R_1,R_2,...,R_p],$ 则： $\tilde{\hat{D}}=[R_p,R_{p-1},...,R_1].$

利用上述符号，L_p+1可表示为如下的迭代表达式：

$L_{p+1}=\left[\begin{array}{cc}{P}_p& {\tilde{a}}_p\\ r_p& L_p\end{array}\right]---\left(12\right)$

上式中：

P_p＝R₀，r_p为由R₁、R₂至R_p排成1列所构成的分块矩阵，根据上面定义的算符可知，

L_p仍为块Toeplitz矩阵，还可继续按照式(12)递归分解，直至L₁＝R₀。

(S5)通过迭代过程求解式(9)的方程，进而得到空时权值：

设则可得到如下的以为未知数的方程：

${\widetilde{\mathrm{\β}}}_p={\widetilde{\mathrm{\α}}}_p{\tilde{L}}_p---\left(13\right)$

易知，当p＝1时，上述方程只涉及对N×N维矩阵的求逆，计算量小，当p＝K时，以上方程即为式(10)的方程。依据文献“H.Akaike.BlockToeplitz Matrix Inversion[J].SIAM J.Appl.Math，1973，24(2)：234-241.”中的理论，可通过一个迭代过程，以较小的计算量求解。在迭代求解得到之后，可得到空时权值为

利用式(7)所示的采样协方差矩阵的共轭对称性，迭代求解空时权值的计算量可进一步缩小，其具体过程如下：

(S51)定义表述迭代过程所需的变量和算符：

设I与Z分别为N×N的单位矩阵与零矩阵。

在第p步迭代中，将方程式(13)的解拆分为p个1×N的向量，即 ${\widetilde{\mathrm{\α}}}_p=[A_1^{\left(p\right)},A_2^{\left(p\right)},...,A_p^{\left(p\right)}].$

对于t＝1，2，…，K-1定义迭代初始化及后续各步骤中的中间变量，具体包括：

1)N×N的矩阵Q(t)、S(t)与G(t)。Q(t)、S(t)的逆矩阵分别表示为Q^-1(t)、S^-1(t)，G(t)的共轭转置表示为G^H(t)；

2)由t个N×N的子矩阵块排成一行得到的矩阵I_s(t)与F_q(t)。对于此两个矩阵序列，定义算符(·)_n表示取其中的第n个子块，例如：(I_s(t))₁表示取I_s(t)的第1个子块，(F_q(t))_t-1表示取F_q(t)中的第t-1个子块。

(S52)迭代初始化：

求解p＝1时，式(13)的方程，具体计算过程如下：

$A_1^{\left(1\right)}=C_1{R}_0^{-1}---\left(14\right)$

求解p＝2时，式(13)的方程，具体计算过程如下：

$I_s\left(1\right)=-R_1{R}_0^{-1},F_q\left(1\right)=-R_{-1}{R}_0^{-1}$

Q^-1(1)＝R₀-F_q(1)I_s(1)R₀，S^-1(1)＝R₀-I_s(1)F_q(1)R₀

                                         (15)

Q(1)＝(Q^-1(1))^-1，S(1)＝(S^-1(1))^-1

$A_2^{\left(2\right)}=(C_2-A_1^{\left(1\right)}{R}_1)Q\left(1\right),A_1^{\left(2\right)}=A_1^{\left(1\right)}+A_2^{\left(2\right)}{F}_q\left(1\right)$

(S53)迭代求解p＝3至p＝K时，式(13)的方程：

依次遍历h＝1，2，…，K-2，执行如下的迭代步骤：

首先，按照下式计算中间变量：

$G\left(h\right)=(I_s\left(h\right){\hat{r}}_h+R_{h+1})$

I_s(h+1)＝[I_s(h)，Z]-G(h)Q(h)[F_q(h)，I]

F_q(h+1)＝[Z，F_q(h)]-G^H(h)S(h)[I，I_s(h)]

Q^-1(h+1)＝Q^-1(h)-(F_q(h+1))₁(I_s(h+1))_h+1Q^-1(h)    (16)

S^-1(h+1)＝S^-1(h)-(I_s(h+1))_h+1(F_q(h+1))₁S^-1(h)

Q(h+1)＝(Q^-1(h+1))^-1，S(h+1)＝(S^-1(h+1))^-1

其次，求解p＝h+2时，方程式(13)的解，具体计算公式如下：

1)计算解向量中的最后一个1×N子块，即计算：

$A_{h+2}^{(h+2)}=(C_{h+2}-\underset{l=1}{\overset{h+1}{\mathrm{\Σ}}}{A}_l^{(h+1)}{R}_{h+2-l})Q(h+1)---\left(17\right)$

2)计算解向量的前h+1个1×N子块，即对于l＝1，2，…，h+1按照下式进行计算：

$A_l^{(h+2)}=A_l^{(h+1)}+A_{h+2}^{(h+2)}{\left(F_q(h+1)\right)}_l---\left(18\right)$

迭代结束后得到的即为式(10)的解。

本发明具备以下有益效果：本发明首先利用在权值更新频度较高时，信号与干扰的平稳性假设，将采样协方差矩阵表示为块Toeplitz矩阵形式，由此首先简化了采用协方差矩阵的计算。然后利用采样协方差矩阵同时为块Toeplitz矩阵和Hermite矩阵的结构特点，本发明提出了迭代求解空时权值的低复杂度方法。与传统的SMI方法相比，在采用相同抗干扰准则的前提下，二者的性能相同，但传统方法的计算复杂度为O[N³ K³]，本发明方法的计算复杂度降低为O[N³ K²]。注意到通常的导航抗干扰天线阵中阵元个数N较小，故求解空时权值的计算量主要取决于时域抽头个数K。本发明方法的计算量与K²成正比，因此随着时域抽头数K的增加，本发明方法相比传统SMI方法减少的计算量越多。较低的计算量一方面可以减少空时滤波器的实现代价；另一方面可加快空时权值的更新频度，以快速适应信号干扰环境的变化。

附图说明

图1表示空时滤波实现结构图；

图2表示本发明中空时权值求解流程图；

图3是有用信号与干扰方向的阵列响应对比图。

具体实施方式

下面，结合附图和具体实施方式，对本发明作进一步说明。

本发明应用于阵列信号处理的空时抗干扰滤波器的滤波权值优化设计中。如附图1所示为具有N个阵元，每个阵元后均有K个时延抽头的空时滤波器结构框图。各天线阵元接收到信号以后经下变频与数模转换以后得到数字中频信号，再通过数字下变频转换为复基带信号。阵元1的复基带信号依次经过K-1个延时单元后得到x₁₁，x₁₂，…，x_1K，阵元2的复基带信号依次经过K-1个延时单元后得到x₂₁，x₂₂，…，x_2K，依此类推阵元N的复基带信号依次经过K-1个延时单元后得到x_N1，x_N2，…，x_NK。与抽头数据x_ab对应的加权系数为各抽头处的数据与加权系数相乘后再全部求和得到阵列的输出y。每经过L个数据采样点，执行附图2所示的步骤S1～S5，以求出空时权值，并将其更新到空时滤波器中。

本发明的目的是减少更新空时权值所需的计算量，只要涉及采样协方差矩阵计算和对其进行求逆运算就可以应用本发明。阵列处理所能获得的性能与具体选择的抗干扰准则有关，下面以采用工程实践中最常用的功率倒置准则为例，用数值仿真来验证本发明的正确性。对于功率倒置准则，c为仅有一个为1，其余全部为0的向量，1对应于参考阵元上的参考抽头。针对导航信号的特点，采用等效阵列增益来衡量抗干扰性能。等效阵列增益描述的是阵列处理对接收机对导航信号解扩后获得的载噪比所产生的影响，取正数表示增强了载噪比，取负数表示对载噪比造成了损耗。

仿真条件设置为：采用4元方阵，阵元间距为半波长，假设信号和干扰都为远场入射，有用信号为GPS的L2频点军用信号，码率为10.23MHz，输入载噪比为43dBHz。数据采样率为(62/3)MHz，采用21000个数据点(对应于约1KHz的权值更新频率)计算采样协方差矩阵(即L＝21000)。设仅有1个有用信号和1干扰信号入射到阵列上，有用信号的方位角为0°，俯仰角为60°，干扰信号的方位角为45°，俯仰角为15°。干扰信号为BPSK调制的宽带干扰，中心频率与信号相同，带宽为20MHz，干信比为60dBc。仿真的时间长度为1400个权值更新周期。

仿真时首先根据仿真条件模拟生成对应于各阵元的模数转换后数据，并将其进一步处理为对应于各阵元的复基带信号。然后通过如下的两步处理来计算等效阵列增益：第一步仅统计出协方差矩阵的估计值，并将得到的协方差矩阵存储下来。第二步根据前一步存储的协方差矩阵，采用本发明提供的步骤S1～S5计算出空时权值，然后再根据空时权值获得空时滤波器的输出，并对此输出信号用软件接收机进行跟踪处理以估计出载噪比，从而得到等效阵列增益。

仿真结果表明，按照本发明方法计算出的空时权值与按照传统SMI方法计算得到的空时权值完全一致。因此与基于功率导致准则的传统SMI方法相比，仅存在采样协方差矩阵计算方式的差异。如表1所示为仿真得到的等效阵列增益与时域抽头个数的关系，从表中可看出，对于不同的时域抽头个数本发明均能有效抑制干扰。

表1等效阵列增益仿真结果

表1中较大的载噪比的损耗主要是由于功率倒置准则未能利用有用信号方向信息导致的。由于仿真时，假设的各阵元及其对应的射频通道完全一致，因此对于不同时域抽头个数，干扰抑制性能基本相同。在时域抽头个数为17时，取一组典型的滤波权值，可计算得到有用信号和干扰信号方法的阵列响应如附图3所示。从附图3可看出，空时率波器在干扰方向的整个带宽范围内产生了35dB以上的零陷，因此有效抑制了干扰信号。由此说明，本发明对采样协方差矩阵的简化计算是正确的，所能获得的抗干扰性能与传统SMI方法相当。

采样协方差矩阵的估计通常由硬件电路以接近实时的速度计算，因此权值更新所能达到的最大频度通常取决于根据采样协方差矩阵计算空时权值的计算量。下面将具体分析本发明步骤S5中计算空时权值的计算量。

分析计算量采用的方法为，先将迭代过程中的计算类型分类，然后再统计每一类运算中所需要的乘法和加法次数，以总的乘加运算次数来表示计算量。从S5的迭代计算过程可知，求解空时权值的计算过程可分解以下五类基本运算：A)N×N矩阵之间的乘法；B)N×N矩阵之间的加减法；C)N×N矩阵的求逆；D)1×N向量与N×N矩阵的乘法；E)1×N向量之间的加减法。通过分析迭代计算过程可得到，在一次空时权值求解过程中，初始化阶段，每一个步迭代过程中所需上述5类基本运算的次数如表2，并据此进一步统计得到了5类运算各自总的执行次数。

表2空时权值求解计算量分析

易知，上述5类基本运算所需的乘法和加法运算总次数分别如下：

A)矩阵乘法的计算量为：

T_mul(N)＝N²*(2N-1)＝O[2N³]    (19)

B)矩阵加减法的计算量可表示为：

T_add(N)＝N²＝O[N²]    (20)

C)由于需要求逆的矩阵均为Hermite矩阵，故利用其对称性可得到其计算量为：

$T_{\mathrm{inv}}\left(N\right)=\frac{2}{3}{N}^3=O\left[\frac{2}{3}{N}^3\right]---\left(21\right)$

D)向量与矩阵乘法的计算量为：

T_vma(N)＝(2N-1)*N＝O[2N²]    (22)

E)向量加减法的计算量为：

T_vpv(N)＝N²＝O[N²]    (23)

综合以上分析可知，本发明单次空时权值计算的总计算量可表示为：

$\begin{array}{c}{T}_{\mathrm{all}}(N,K)=(\frac{3}{2}{K}^2+\frac{3}{2}K-3)T_{\mathrm{mul}}\left(N\right)+\\ (\frac{3}{2}{K}^2-\frac{5}{2}K+1)T_{\mathrm{add}}\left(N\right)\\ +(2K-1)T_{\mathrm{inv}}\left(N\right)\\ +(K^2-K+2)T_{\mathrm{vma}}\left(N\right)+\\ (K^2-K)T_{\mathrm{vpv}}\left(N\right)=O\left[{3N}^3{K}^3\right]\end{array}---\left(24\right)$

上式表明，本发明步骤S5的计算量小于传统的SMI方法，且由于计算量与K²成正比，因此对于越大的时域抽头数K，本发明相对于传统SMI方法减少的计算量越多。

综合数值仿真和计算量分析可知，本发明的空时抗干扰方法在保持与传统SMI类方法具有相同性能的前提下，大幅减少了计算量。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims (2)

1.一种块Toeplitz矩阵低复杂度求逆的空时抗干扰方法，其特征在于，包括以下步骤：

(S1)初始化，将空时权值求解问题转化为线性方程组求解问题；

设第a个阵元后第b个时域抽头处的复基带信号为采样值序列x_ab[j]，其中，j为整数，表示采样时刻；空时滤波器中与之相应的加权值的复共轭为w_ab，将同一采样时刻各抽头处的采样值先按照阵元先后顺序排列，再按照抽头先后顺序排列，构成如下的随机向量：

X＝[x₁₁，x₂₁，…，x_N1，x₁₂，x₂₂，…，x_N2，…，x_1K，x_2K，…，x_NK]^T     (1)

其中，N表示阵元个数，K表示抽头个数；

相应地，将所有的空时权值的复共轭写为如下的向量形式：

W＝[w₁₁，w₁₂，…，w_1K，w₂₁，w₂₂，…，w_2K，…，w_N1，w_N2，…，w_NK]^T      (2)

空时滤波器的输出表示为：

y＝W^HX(3)

其中，T表示矩阵转置，(·)^H表示对矩阵或向量进行共轭转置操作；

阵列各抽头处数据的协方差矩阵为：E(XX^H)，E(·)表示求数学期望运算，将数学期望用求平均值取代，得到作为协方差矩阵的采样协方差矩阵R_xx；

设在一个权值更新周期内有L个采样点，则R_xx的第u行第υ列元素表示如下：

${\left(R_{\mathrm{xx}}\right)}_{\mathrm{u\υ}}=\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{\left(X\right)}_u{\left(X\right)}_{\mathrm{\υ}}^{*}---\left(4\right)$

上式中，(X)_u表示向量X中的第u个元素，(X)_v表示向量X中的第υ个元素，(·)^*表示求复数的转置；

将空时权值转化为如下形式的线性方程组的解：

R_xxW＝c          (5)

上式中，c为由抗干扰优化准则所决定的约束向量；

(S2)将采样协方差矩阵表示为块Toeplitz矩阵的形式，并简化其计算量；

将采样协方差矩阵采用如下的块Toeplitz矩阵形式来计算：

$R_{\mathrm{xx}}=\left[\begin{array}{ccccc}{R}_0& R_1& R_2& ...& R_{K-1}\\ R_{-1}& R_0& R_1& ...& R_{K-2}\\ R_{-2}& R_{-1}& R_0& ...& R_{K-3}\\ \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·\\ R_{-K+1}& R_{-K+2}& R_{-K+3}& ...& R_0\end{array}\right]---\left(6\right)$

其中，R_i，i＝-K+1，…，0，…，K-1为2K-1个N×N方阵；

根据各子块之间满足如下的共轭对称性：

$R_i^H=R_{-i},i=-K+1,...,0,...,K-1$

将计算采样协方差矩阵简化为计算K个N×N的子阵R_m，m＝0，…，K-1；R_m的第u行第υ列元素按照下式计算：

${\left(R_m\right)}_{\mathrm{u\υ}}=\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{x}_{u1}\left[j\right]x_{\mathrm{\υ}(m+1)}^{*}\left[j\right]---\left(8\right)$

(S3)将步骤(S1)中的线性方程组进行简化处理为迭代求解问题；

将空时权值的转置及约束向量的转置均表示为K个1×N的行向量构成的分块矩阵，即：

W^H＝[A₁，A₂，…，A_K] c^H＝[C₁，C₂，…，C_K]     (9)

利用R_xx的共轭对称性，将式(5)的空时权值求解问题表示为：

[C₁，C₂，…，C_K]＝[A₁，A₂，…，A_K]R_xx        (10)

(S4)定义描述迭代求解方程的变量及算符；

设迭代步数为p＝1，2，…，K，L_p+1为第一行子块为R₀、R_-1至R_-p，第一列子块为R₀、R₁至R_p所构成的块Toeplitz矩阵，即：

$L_{p+1}=\left[\begin{array}{ccccc}{R}_0& R_{-1}& R_{-2}& ...& R_{-p}\\ R_1& R_0& R_{-1}& ...& R_{-p+1}\\ R_2& R_1& R_0& ...& R_{-p+2}\\ \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·\\ R_p& R_{p-1}& R_{p-2}& ...& R_0\end{array}\right]---\left(11\right)$

对于构成元素均为N_×N矩阵的任意分块矩阵D，引入两个算符：表示对矩阵进行以子块为单位的转置操作，即中第(u，v)子块，正好是D中的第(υ,u)子块；表示交换矩阵中同一行中各子块的先后顺序，即若设： $\tilde{D}=[R_1,R_2,...,R_p],$ 则： $\stackrel{\widetilde{^}}{D}=[R_p,R_{p-1},...,R_1];$

利用上述符号，L_p-1表示为如下的迭代表达式：

$L_{p+1}=\left[\begin{array}{cc}{P}_p& {\tilde{a}}_p\\ r_p& L_p\end{array}\right]---\left(12\right)$

其中，P_p＝R₀，r_p为由R₁、R₂至R_p排成1列所构成的分块矩阵，根据上面定义的算符得到， ${\tilde{r}}_p=[R_1,R_2,...,R_p];$

L_p仍为块Toeplitz矩阵，继续按照式(12)递归分解，直至L₁＝R₀；

(S5)通过迭代过程求解得到空时权值；

设则得到如下的以为未知数的方程：

${\widetilde{\mathrm{\β}}}_p={\widetilde{\mathrm{\α}}}_p{\tilde{L}}_p---\left(13\right)$

当p＝1时，方程(13)为对N×N维矩阵的求逆，当p＝K时，方程(13)为式(10)的方程，通过迭代过程求解，根据迭代求解得到结果后，计算得到空时权值为

2.如权利要求1所述的一种块Toeplitz矩阵低复杂度求逆的空时抗干扰方法，其特征在于，所述步骤(S5)中的迭代过程为：

(S51)定义表述迭代过程所需的变量和算符：

设I与Z分别为N×N的单位矩阵与零矩阵；

在第p步迭代中，将方程式(13)的解拆分为p个1×N的向量，即 ${\widetilde{\mathrm{\α}}}_p=[A_1^{\left(p\right)},A_2^{\left(p\right)},...,A_p^{\left(p\right)}];$

对于t＝1，2，…，K-1，定义迭代初始化及运算步骤中的中间变量，具体包括：

1)N×N的矩阵Q(t)、S(t)与G(t)；Q(t)、S(t)的逆矩阵分别表示为Q^-1(t)、S^-1(t)；G(t)的共轭转置表示为G^H(t)；

2)由t个N×N的子矩阵块排成一行得到的矩阵I_s(t)与F_q(t)；对于此两个矩阵序列，定义算符(·)_n表示取其中的第n个子块，如：(I_s(t))₁表示取I_s(t)的第1个子块，(F_q(t))_t-1表示取F_q(t)中的第t-1个子块；

(S52)迭代初始化；

求解p＝1时，式(13)的方程，具体计算过程如下：

$A_1^{\left(1\right)}=C_1{R}_0^{-1}---\left(14\right)$

求解p＝2时，式(13)的方程，具体计算过程如下：

$\begin{array}{c}{I}_s\left(1\right)=-R_1{R}_0^{-1},F_q\left(1\right)=-R_{-1}{R}_0^{-1}\\ Q^{-1}\left(1\right)=R_0-F_q\left(1\right)I_s\left(1\right)R_0,S^{-1}\left(1\right)=R_0-I_s\left(1\right)F_q{R}_0\\ Q\left(1\right)={\left(Q^{-1}\left(1\right)\right)}^{-1},S\left(1\right)={\left(S^{-1}\left(1\right)\right)}^{-1}\\ A_2^{\left(2\right)}=(C_2-A_1^{\left(1\right)}{R}_1)Q\left(1\right),A_1^{\left(2\right)}=A_1^{\left(1\right)}+A_2^{\left(2\right)}{F}_q\left(1\right)\end{array}---\left(15\right)$

(S53)迭代求解p＝3至p＝K时，式(13)的方程：

依次遍历h＝1，2，…，K-2，执行如下的迭代步骤：

首先，按照下式计算中间变量：

$\begin{array}{c}G\left(h\right)=(I_s\left(h\right){\hat{r}}_h+R_{h+1})\\ I_s(h+1)=[I_s\left(h\right),Z]-G\left(h\right)Q\left(h\right)[F_q\left(h\right),I]\\ F_q(h+1)=[Z,F_q\left(h\right)]-G^H\left(h\right)S\left(h\right)[I,I_s\left(h\right)]\\ Q^{-1}(h+1)=Q^{-1}\left(h\right)-{\left(F_q(h+1)\right)}_1{\left(I_s(h+1)\right)}_{h+1}{Q}^{-1}\left(h\right)\\ S^{-1}(h+1)=S^{-1}\left(h\right)-{\left(I_s(h+1)\right)}_{h+1}{\left(F_q(h+1)\right)}_1{S}^{-1}\left(h\right)\\ Q(h+1)={\left(Q^{-1}(h+1)\right)}^{-1},S(h+1)={\left(S^{-1}(h+1)\right)}^{-1}\end{array}---\left(16\right)$

其次，求解p＝h+2时，方程式(13)的解，具体计算公式如下：

1)计算解向量中的最后一个1×N子块，即计算：

$A_{h+2}^{(h+2)}=(C_{h+2}-\underset{l=1}{\overset{h+1}{\mathrm{\Σ}}}{A}_l^{(h+1)}{R}_{h+2-l})Q(h+1)---\left(17\right)$

2)计算解向量的前h+1个1×N子块，即对于l＝1，2，…，h+1按照下式进行计算：

$A_l^{(h+2)}=A_l^{(h+1)}+A_{h+2}^{(h+2)}{\left(F_q(h+1)\right)}_l---\left(18\right)$

迭代结束后，得到的]为式(10)的解。