CN116915552B - Irs辅助级联mimo双稀疏信道状态信息估计方法 - Google Patents

Abstract

本发明涉及IRS辅助级联MIMO双稀疏信道状态信息估计方法，方法包括以下步骤：S1、获取经角度过采样后的级联信道的表达式，再将该表达式向量化；S2、建立状态空间模型；S3、初始化滤波参数；S4、获取当前导频符号的观测向量；S5、由卡尔曼滤波方法计算一步预测和估计角度增益值；S6、初始化带有范数约束的卡尔曼滤波的参数；S7、用带有范数约束的卡尔曼滤波方法一步预测和估计角度增益值；S8、判断两次相邻的估计值数的差值绝对值是否小于阈值，进行迭代计算；S9、更新符号为当前符号的下一个符号，重复S4～S8，滤波达到稳态后，得到瞬时的级联信道。与现有技术相比，本发明具有适用性强，无需先验信息等优点。

Description

IRS辅助级联MIMO双稀疏信道状态信息估计方法

技术领域

本发明涉及无线信道状态信息估计的技术领域，尤其是涉及一种IRS辅助级联MIMO双稀疏信道状态信息估计方法。

背景技术

智能反射表面(IRS)是带有大量被动反射元件的超表面。它可以改变入射信号的相位，起到调节无线传输环境和提高频谱效率的作用，在直射链路受到遮挡时增强无线覆盖。在引入IRS的无线传输系统中获取辅助级联的信道状态信息(CSI)是性能提升中必不可少的一环。

由于毫米波波段和太赫兹波段的无线信道通常只有少数的弥散路径，IRS辅助级联信道的每跳信道都是角度域和延迟域的双稀疏信道，采用压缩感知理论框架下的稀疏信号表示和恢复来估计毫米波和太赫兹波段级联信道已获得一定的研究。

考虑到级联信道的高维性和稀疏性，结合平行因子分解模型，学者们已提出一系列交替的稀疏信道估计算法，所使用的压缩感知类(Compressive Sensing，CS)算法包括两大类。第一类为在稀疏度已知时，以正交匹配追踪(OMP)算法为代表的稀疏重构算法，包括双交替OMP(BAOMP)算法和双交替正则化OMP(BAROMP)算法；第二类为稀疏度未知时，通过设置步长和适当的终止条件估计稀疏度和填充支集，给出了以稀疏度自适应匹配追踪(SAMP)算法为代表的稀疏重构算法，这些贪婪类稀疏恢复算法通常仅利用了导频音上的单测量(观测)向量，稀疏恢复能力有限，在级联信道获取中的性能还有待提高。再有，基于张量分解的交替稀疏恢复方法尚存在着弊端，由于张量分解的不唯一性，在分解后不能保证两段信道的稀疏位置与原信道的稀疏位置一样，限制了要求稀疏度已知的恢复算法的应用。

综上所述，现有的稀疏信道估计算法在字典矩阵高相关和信道秩亏严重时，恢复性能差且需要稀疏度先验信息，适用性不强。

发明内容

本发明的目的就是为了克服需要先验信息导致适用性不强的问题而提供的一种IRS辅助级联MIMO双稀疏信道状态信息估计方法。

本发明的目的可以通过以下技术方案来实现：

一种IRS辅助级联MIMO双稀疏信道状态信息估计方法，方法包括以下步骤：

S1、对于配置N_T根发送天线的通信发送端、带N_I个反射元件的IRS和配置N_R根接收天线的通信接收端所组成的级联传输系统，为通信发送端与IRS端之间的信道H₁，构造接收角度域字典矩阵D_1R和发送角度域字典矩阵D_1T，同时为通信接收端与IRS端之间的信道H₂，构造接收角度域字典矩阵D_2R和发送角度域字典矩阵D_2T，基于上述字典矩阵，得到经角度过采样后的级联信道的表达式，再将该表达式向量化，向量化的表达式中包括待求的级联角度增益向量h；

S2、对S1的级联角度增益向量建立任意的第n个符号对应的状态空间模型h(n)；

S3、初始化级联角度增益向量h为h(0)，初始化级联角度增益向量h的估计误差方差为P(0)，设置状态空间模型中更新噪声的协方差矩阵Q，根据信噪比设置观测噪声方差R；

S4、通过第s个符号的带有保护间隔迫零块对角导频X(s)获得导频音观测矩阵Y(s)，并进行向量化，得到当前观测向量y(s)，所述观测向量y(s)和符号序号s一一对应；

S5、对于当前观测向量，由卡尔曼滤波方法计算角度增益向量的一步预测h(s|s-1)，根据协方差矩阵Q计算一步预测估计方差P(s|s-1)，根据观测噪声方差R计算滤波增益K(s)，根据滤波增益和角度增益向量的一步预测得到当前估计值h(s)，根据滤波增益和一步预测估计方差得到当前估计误差方差P(s)；

S5中第一次采用卡尔曼滤波方法计算时的初始值采用S3的初始化的h(0)和P(0)，其他次数的初始值采用S8的当前第s个符号的级联角度增益向量和估计方差；

S6、构造l₁范数的伪测量方程，设置l₁范数的误差方差，把当前估计值h(s)作为l₁范数迭代修正的初始值h_l(0)，把当前的估计误差方差P(s)作为迭代初始估计误差方差P_l(0)；

S7、用卡尔曼滤波方法计算一步预测h_l(τ|τ-1)，计算一步预测估计误差方差P_l(τ|τ-1)，然后计算迭代滤波增益k_l(τ)，计算当前估计值h_l(τ)和当前估计误差方差P_l(τ)，其中τ为迭代次数；

S7中第一次采用卡尔曼滤波方法计算时的初始值采用S6的初始值h_l(0)和初始估计误差方差P_l(0)；

S8、判断两次相邻的估计值h_l(τ)和h_l(τ-1)的l₁范数的差值绝对值是否小于阈值，若是则停止迭代，反之则判断是否达到最大迭代次数，若是则停止迭代，反之继续执行S7进行迭代；

S8中，停止迭代后，把已进行l₁范数迭代修正的h_l(τ)作为当前第s个符号的级联角度增益向量，把停止迭代时的P_l(τ)作为当前第s个符号的估计误差方差；

S9、更新符号为当前符号s的下一个符号，重复S4～S8，滤波达到稳态后，得到稀疏的角度增益向量中的角度支集和相应的角度增益，根据角度域字典得到瞬时的矩阵向量化的信道，通过矩阵化得到瞬时的级联信道H(s)。

进一步地，S1的向量化的表达式为：

其中，H表示经角度过采样后的级联信道，g(s)表示IRS在第s个符号上的相移向量，A_R，I和A_I，T分别为在字典导向矢量方向上通信接收端和IRS端之间的角度增益矩阵以及IRS端和通信发送端之间的角度增益矩阵，符号和⊙分别表示克罗内克积和Khatri-Rao积，符号(·)*，(·)^T和(·)^H分别表示共轭、转置和共轭转置，为待求的级联角度增益向量。

进一步地，S5的角度增益向量的一步预测h(s|s-1)和一步预测估计方差P(s|s-1)为：

h(s|s-1)＝ψh(s-1)

P(s|s-1)＝ψP(s-1)+(1-ψ²)Q

其中，ψ为状态空间模型中由最大多普勒频移f_D确定的参数，P(s-1)为上一个符号对应的估计误差方差，Q为更新噪声的协方差矩阵。

进一步地，滤波增益K(s)为：

K(s)＝P(s|s-1)A^H(s)(A(s)P(s|s-1)A^H(s)+R)^-1

其中，R为观测噪声方差，X(s)为带有保护间隔迫零块对角导频符号，(·)^*，(·)^T和(·)^H分别表示共轭、转置和共轭转置；

当前估计值h(s)和当前估计误差方差P(s)为：

h(s)＝h(s|s-1)+K(s)(y(s)-A(s)h(s|s-1))

P(s)＝P(s|s-1)-K(s)A(s)P(s|s-1)

其中，h(s|s-1)为角度增益向量的一步预测，P(s|s-1)为一步预测估计误差方差。

进一步地，S6的伪测量方程为：

0＝||h(s)||₁-η(τ)∈

其中，η为自适应因子，∈为一个方差固定的随机数，τ为迭代次数。

进一步地，自适应因子为：

η(τ)＝1-η(0)exp(-μτ/T)

其中，η(0)和μ是形状参数，T为一个固定的正常数，τ为迭代次数。

进一步地，S7的一步预测η_l(τ|τ-1)和一步预测估计误差方差P_l(τ|τ-1)为：

h_l(τ|τ-1)＝h_l(τ-1)

P_l(τ|τ-1)＝P_l(τ-1)

其中，h_l(τ-1)为上一迭代次数时的估计值，P_l(τ-1)为上一迭代次数时的估计误差方差；

迭代滤波增益k_l(τ)为：

k_l(τ)＝P_l(τ|τ-1)c^H(τ)(c(τ)P_l(τ|τ-1)c^H(τ)+r(τ))^-1

其中，c(τ)为由次梯度向量构成的雅可比矩阵，r(τ)＝η²(τ)r_∈为l₁范数的测量方差，r_∈为随机数∈的固定方差，η为自适应因子。

进一步地，当前估计值h_l(τ)和当前估计误差方差P_l(τ)为：

h_l(τ)＝h_l(τ|τ-1)+k_l(τ)(0-c(τ)h_l(τ|τ-1))

P_l(τ)＝P_l(τ|τ-1)-k_l(τ)c(τ)P_l(τ|τ-1)

其中，k_l(τ)为迭代滤波增益。

进一步地，滤波达到稳态的判断条件为：

基于相邻两个符号的S8得到的估计误差方差的主对角线元素和的绝对差是否小于阈值进行判断，若是，则达到稳态，反之未达到稳态。

进一步地，接收角度域字典矩阵D_1R和发送角度域字典矩阵D_1T分别为：

接收角度域字典矩阵D_2R和发送角度域字典矩阵D_2T分别为：

其中，a_1r，j和a_1t，j，j＝0，...，β₁-1是信道H₁在角度网格点上的导向矢量，a_2r，j和a_2t，j，j＝0，...，β₂-1是信道H₂在角度网格点上的导向矢量。

与现有技术相比，本发明具有以下有益效果：

本发明利用角度域过采样字典和带有保护间隔迫零块对角导频，通过序列的导频音观测和稀疏角度增益向量的时域相关性，用带有范数约束的卡尔曼滤波方法重构静态或动态的稀疏角度增益向量和相应的角度支集和角度增益向量，可用较少的导频开销得到准确的信道状态信息。本发明对角度域过采样字典没有限制，可在静态信道场景和动态信道场景下获取角度支集，增强了适用性，同时不需要任何关于稀疏角度增益向量的稀疏度先验信息，对静态信道的角度增益和支集估计精度高，对动态信道可估计出角度支集并跟踪角度增益。

附图说明

图1为本发明的流程图；

图2为本发明的应用场景示意图，其中图2(a)为固定位置的通信发送端(如用户)到IRS再到固定位置的通信接收端(如基站)的传输场景，图2(b)为移动的通信发送端(如车辆)到IRS再到固定的通信接收端(如基站)的应用场景；

图3为本发明的角度域角度增益矩阵的稀疏结构示意图；

图4为本发明的带有保护间隔迫零块对角导频图案示意图；

图5为本发明在静态场景1(N_I＝4，4×4维MIMO，L₁＝2，L₂＝3)下，角度增益的归一化均方误差(NMSE)随信噪比(SNR)的变化曲线；

图6为本发明在静态场景1(N_I＝4，4×4维MIMO，L₁＝2，L₂＝3)下，对信道复衰落幅值的跟踪曲线；

图7为本发明在静态场景1(N_I＝4，4×4维MIMO，L₁＝2，L₂＝3)下，对角度支集和角度增益的估计情况；

图8为本发明在静态场景1(N_I＝4，4×4维MIMO，L₁＝2，L₂＝3)下，当SNR＝10dB时归一化均方误差(NMSE)随导频音观测向量个数的变化曲线；

图9为本发明在动态场景1(N_I＝4，4×4维MIMO，L₁＝2，L₂＝3)下，当SNR＝10dB和导频音观测向量个数为20时，归一化均方误差(NMSE)随ψ的变化曲线；

图10为本发明在动态场景1(N_I＝4，4×4维MIMO，L₁＝2，L₂＝3)下，在SNR＝10dB、导频音观测向量数为20和ψ＝0.9时，对一个非零的角度增益h₁的跟踪曲线；

图11为本发明在动态场景1(N_I＝4，4×4维MIMO，L₁＝2，L₂＝3)下，在SNR＝10dB和ψ＝0.9时，所估计的角度增益向量的l₁范数随导频音观测向量个数的变化曲线；

图12为本发明在静态场景2(N_I＝4，4×4维MIMO，L₁＝1，L₂＝1)下，角度增益的归一化均方误差(NMSE)随SNR的变化曲线；

图13为本发明在静态场景2(N_I＝4，4×4维MIMO，L₁＝1，L₂＝1)下，对信道复衰落幅值的跟踪曲线；

图14为本发明在静态场景2(N_I＝4，4×4维MIMO，L₁＝1，L₂＝1)下，对角度支集和角度增益的估计情况；

图15为本发明在静态场景2(N_I＝4，4×4维MIMO，L₁＝1，L₂＝1)下，当SNR＝10dB时归一化均方误差(NMSE)随导频音观测向量个数的变化曲线；

图16为本发明在动态场景2(N_I＝4，4×4维MIMO，L₁＝1，L₂＝1)下，当SNR＝10dB和导频音观测向量个数为20时，归一化均方误差(NMSE)随ψ的变化曲线；

图17为本发明在动态场景2(N_I＝4，4×4维MIMO，L₁＝1，L₂＝1)下，在SNR＝10dB、导频音观测向量数为20和ψ＝0.9时，对一个非零的角度增益h₁的跟踪曲线；

图18为本发明在动态场景2(N_I＝4，4×4维MIMO，L₁＝1，L₂＝1)下，在SNR＝10dB和ψ＝0.9时，所估计的角度增益向量的l₁范数随导频音观测向量个数的变化曲线；

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

名词解释：

智能反射表面——IRS；

多输入多输出——MIMO。

本发明提出一种IRS辅助级联MIMO双稀疏信道状态信息估计方法，方法的流程图如图1所示。本发明的应用场景包括如图2所示的场景，它包括(a)中所示的固定位置的通信发送端(如用户)到IRS再到固定位置的通信接收端(如基站)的传输场景，也包括(b)中所示的移动的通信发送端(如车辆)到IRS再到固定的通信接收端(如基站)的应用场景，还可包括与(a)和(b)所示类似的下行链路应用场景。在(a)和(b)中，通信发送端的多天线和IRS元件间的MIMO信道H₁(s)可为时不变的静态信道或时变的动态信道，从IRS元件到通信接收端天线的MIMO信道H₂(s)为时不变信道。

本发明的IRS辅助级联MIMO双稀疏信道模型指毫米波或太赫兹波段的方向阵列信道模型，其中，通信发送端和通信接收端的天线阵列既可为平面阵，也可为线阵，IRS元件阵列一般为平面阵，在特殊情况下可为线阵。不失一般性。设两跳信道的每一跳均为毫米波平面阵列天线信道，信道的角度域参数与路径数、路径波束的离开角、到达角以及相应的导向矢量有关，可表示为：

在式(1)和式(2)中，L₁和L₂表示两段信道的路径数，和表示路径角度增益，δ(·)为脉冲成形滤波器响应，和分别为路径延迟，和以及 分别表示信号入射(上标r)及出射(上标t)的方位角和仰角，导向矢量a_r和a_t可统一写为：

式中，0≤m<N_y，0≤n<N_z是平面阵元在y轴和z轴的索引序号，k＝2π/λ，λ是波长，d是天线间隔，d＝0.5λ，(·)^T表示转置。

当角度参数和以及和与相应的角度增益和乃至路径数L₁和L₂已知，就可以确定一个角度域稀疏的毫米波信道模型，所说的延迟域稀疏，是指多路径响应在延迟域上是可分离的。为分离每个路径响应拍，可利用带有保护间隔迫零块对角导频图案，通过插入迫零保护间隔来分离多拍响应为多个单拍响应的和，从而避免了符号间干扰。毫米波或太赫兹波信道双稀疏特性里的延迟域稀疏特性是隐藏的，可利用这个隐藏的稀疏性来设计导频图案。

在波束方向未知即依赖于离开角和到达角的导向矢量未知时，需要用空间角度网格来划分空间角度并探测响应。对通信发送端和IRS端之间的信道H₁用β₁个角度网格划分，每个网格对应一对发送端的方位角和仰角以及一对接收端的方位角和仰角，再用式(3)可分别得到在每个网格上的导向矢量，然后得到信道H₁在角度域过采样的字典矩阵，包括IRS为接收端时的字典矩阵和通信发送端的字典矩阵同样地，对通信接收端和IRS端间的信道H₂用β₂个角度网格划分，每个网格也对应一对发送端的方位角和仰角以及一对接收端的方位角和仰角，得到H₂在角度域过采样的字典矩阵，包括通信接收端的字典矩阵和IRS为发送端时的字典矩阵

设对应这些角度域字典的角度增益分别为稀疏矩阵和通过角度域字典，级联信道可以写为：

式中，g(s)为在符号s上的N_I×1维IRS相移向量，diag(·)表示向量对角化算子，A_R，I和A_I，T分别为字典矩阵所对应的通信接收端与IRS元件间的稀疏角度增益矩阵和IRS元件和通信发送端间的稀疏角度增益矩阵，(·)^H表示共轭转置。

根据克罗内克(Kronecker)积(用符号表示)和Khatri-Rao积(用符号⊙表示)的性质，式(4)可以写为：

式中，vec(·)表示按列将矩阵向量化的算子，(·)*表示共轭操作；为待求的稀疏角度增益向量，并记为

在考虑通信发送端或通信接收端的移动性时，h可建模为状态空间模型：

式中，ψ为角度增益向量的时域相关系数，它由Jakes信道模型ψ＝J₀(2πf_DT_s)确定，其中，J₀(·)表示第一类零阶贝塞尔函数，f_D为通信中的移动一方相对于IRS的最大多普勒频移，T_s为符号周期，当通信发送端和通信接收端均相对于IRS静止，ψ＝1，角度增益向量的状态空间模型退化为时不变信道模型h(s)＝h(s-1)，v(s-1)为更新噪声。

在这个状态空间模型中，h中的零元素不随符号变化总是保持为零，当ψ＝1时为静态的稀疏信道模型，当ψ<1时为动态的稀疏信道模型。

本发明通过角度域过采样字典和带有保护间隔迫零块对角导频，在卡尔曼滤波框架下用序列的导频音观测重构静态或动态的稀疏角度增益向量，利用结构化角度增益向量的时域相关性，用时序递推方式准确估计出角度支集和角度增益向量，可用较少的导频开销得到准确的信道状态信息。本发明的主要思路是，为了增强适用场景和减少导频开销，同时减小角度域信道状态信息估计的计算量，对MIMO级联信道的两段信道采用不加限制的角度网格字典并构造级联后的角度字典矩阵，在卡尔曼滤波框架下，利用导频音序贯测量值，用带有l₁范数约束的卡尔曼滤波器估计角度支集和角度增益。它的优点有：它对角度字典没有任何限制，可在静态信道场景和动态信道场景下获取角度域的支集，因而增强了本发明的适用性；它不需要任何关于稀疏度的先验信息，可用少量的导频音观测精确地获取结构化稀疏信道，在静态环境下对角度增益的估计精度高，在动态环境下可获取较精确的角度支集。

本发明的方法包括以下步骤：

A、对配置N_T根发送天线的通信发送端、带N_I个反射元件的IRS和配置N_R根接收天线的通信接收端所组成的级联传输系统，为通信发送端与IRS端之间的信道H₁，构造接收角度域字典矩阵和发送角度域字典矩阵为通信接收端与IRS端之间的信道H₂，构造接收角度域字典矩阵和发送角度域字典矩阵其中β₁和β₂分别是在角度域上划分H₁和H₂两段信道的角度网格总数，其中，a_1r，j和a_1t，j，j＝0，...，β₁-1是信道H₁在角度网格点上的导向矢量，a_2r，j和a_2t，j，j＝0，...，β₂-1是信道H₂在角度网格点上的导向矢量；由IRS在第s个符号的相移向量g(s)，对通信发送端与通信接收端之间的N_R×N_T维级联信道H＝H₂diag(g(s))H₁，经角度过采样后表示为再将该矩阵形式的信道向量化为其中，下标T、I、R分别表示通信发送端、IRS端和通信接收端，vec(·)表示矩阵的向量化算子，diag(·)表示向量的对角化算子，A_R，I和A_I，T分别为在字典导向矢量方向上通信接收端和IRS端之间的角度增益矩阵以及IRS端和通信发送端之间的角度增益矩阵，符号和⊙分别表示Kronecker积和Khatri-Rao积，符号(·)^*，(·)^T和(·)^H分别表示共轭、转置和共轭转置，为维待求的级联角度增益向量。

步骤A的通信接收端和通信发送端，可为位置固定的通信两端，或一端移动、一端固定的通信两端，对MIMO收发天线数量和IRS元件数量没有限制，过采样的角度网格数可设置为2～4倍的收发阵列天线数，当收发天线数较大时，可设置为1倍的收发阵列天线数，即β₁＝N_T和β₂＝N_R，对平面阵和线阵MIMO天线，角度网格划分可依据角度分辨率均匀采样，也可以随机采样；若在某个角度字典上未得到幅值较大的角度增益，可在粗网格上再细分角度网格来得到角度域上的信道状态信息。

B、为角度增益向量h建立状态空间模型其中，s表示符号索引，ψ由最大多普勒频移f_D确定，v(s-1)为更新噪声，h中的零元素不随符号变化始终保持为零。

步骤B的角度增益向量的状态空间模型中的ψ，由Jakes信道模型ψ＝J₀(2πf_DT_s)确定，其中J₀(·)表示第一类零阶贝塞尔函数，f_D为通信中的移动一方相对于IRS的最大多普勒频移，T_s为符号周期，当通信发送端和通信接收端均相对于IRS静止，ψ＝1，角度增益向量的状态空间模型退化为时不变信道模型h(s)＝h(s-1)。

C、初始化角度增益向量h为h(0)，初始化h的估计误差方差为P(0)，设置更新噪声v(s-1)的协方差矩阵Q，根据信噪比设置观测噪声方差R。

D、通过第s个符号的带有保护间隔迫零块对角导频X(s)，获得N_R×N_T维导频音观测矩阵Y(s)＝H(s)X(s)+W(s)，其中，W(s)为观测噪声矩阵，并得到向量化的观测方程y(s)＝A(s)h(s)+w(s)，其中，y(s)＝vec(Y(s))为N_RN_T×1维观测向量，w(s)＝vec(W(s))为N_RN_T×1维观测噪声向量。

E、由卡尔曼滤波方法计算角度增益向量的一步预测为h(s|s-1)＝ψh(s-1)，计算一步预测估计方差为P(s|s-1)＝ψP(s-1)+(1-ψ²)Q，计算增益为K(s)＝P(s|s-1)A^H(s)(A(s)P(s|s-1)A^H(s)+R)^-1，计算当前估计值为h(s)＝h(s|s-1)+K(s)(y(s)-A(s)h(s|s-1))，计算当前估计误差方差为P(s)＝P(s|s-1)-K(s)A(s)P(s|s-1)。

F、构造l₁范数的伪测量方程为0＝||h(s)||₁-η(τ)∈，其中，η(τ)为一个随迭代次数τ单调递增的自适应因子，当τ趋于无穷时，η(τ)趋于1，∈为一个方差固定的随机数，记∈的方差为r_∈，将l₁范数的误差方差设置为r＝η²(τ)r_∈，把当前估计值h(s)作为l₁范数迭代修正的初始值h_l(0)，把当前的估计误差方差P(s)作为初始估计误差方差P_l(0)。

步骤F的l₁范数伪测量方程，为方差自适应的伪测量方程，它引入了自适应因子η(τ)来加快l₁范数收敛，并且收敛过程可用一个指数增长的η(τ)＝1-η(0)exp(-μτ/T)来调整，其中，η(0)和μ是形状参数，T为一个固定的正常数。

G、用卡尔曼滤波方法计算一步预测为h_l(τ|τ-1)＝h_l(τ-1)，计算一步预测估计误差方差为P_l(τ|τ-1)＝P_l(τ-1)，在当前估计值h_l(τ)上计算l₁范数的次梯度：

得到的次梯度向量c(τ)，它是将非线性l₁范数伪测量方程线性化的雅可比矩阵，然后计算迭代滤波增益为k_l(τ)＝P_l(τ|τ-1)c^H(τ)(c(τ)P_l(τ|τ-1)c^H(τ)+r(τ))^-1，计算当前估计值为h_l(τ)＝h_l(τ|τ-1)+k_l(τ)(0-c(τ)h_l(τ|τ-1))，计算当前估计误差方差为P_l(τ)＝P_l(τ|τ-1)-k_l(τ)c(τ)P_l(τ|τ-1)。

步骤G的l₁范数的次梯度是复数角度增益向量的l₁范数的次梯度，当角度增益为实数时，l₁范数对非零元素的次梯度为符号函数，对零元素的次梯度为(-1，1)的随机值；所述的步骤G的一步预测h_l(τ|τ-1)＝h_l(τ-1)和一步预测估计误差方差P_l(τ|τ-1)＝P_l(τ)仅为赋值操作，不含有计算量。

H、判断相邻两次估计值的l₁范数的差值绝对值是否小于一个小的正数，即判断是否满足|||h_l(τ)||1-||h_l(τ-1)||₁|<ζ，其中，|·|表示绝对值，ζ为一个小的正数，若满足停止迭代，若不满足再判断是否达到最大迭代次数τ_max，若满足停止迭代，若不满足则进入步骤G继续迭代；停止迭代后，把已进行l₁范数迭代修正的h_l(τ)作为当前h(s)，把P_l(τ)作为当前的P(s)。

步骤H的相邻两次估计值的l₁范数的差值绝对值小于一个小的正数条件是迭代的一个终止条件，为保证算法可靠进行，可设置达到最大迭代次数τ_max为另一个终止条件，两个终止条件的其中一个条件满足时，l₁范数修正的迭代过程停止。

I、更新符号为当前符号s的下一个符号，重复步骤D、E、F、G、H进行估计和跟踪h(s)，滤波达到稳态后，得到稀疏的角度增益向量h(s)中的角度支集即非零元素在向量中的位置和相应复角度增益，由角度字典得到瞬时的 再通过矩阵化得到瞬时的级联信道H(s)。

步骤I的滤波达到稳态，可用相邻两个符号的估计误差方差P(s)和P(s-1)的主对角线元素和的绝对差是否已是小量来判定。

下面以N_R＝N_I＝N_T＝4和L₁＝2及L₂＝3时，如附图3所示为例，对本发明的方法进行说明。

在N_R＝N_I＝N_T＝4和L₁＝2及L₂＝3时，一种采用1倍阵列尺寸的角度字典所对应的角度增益矩阵的稀疏结构，它的任一行或任一列要么全是0，要么只有一个非零，该稀疏结构也可为其他结构，本发明方法对稀疏结构没有限制。

如附图4所示为带有保护间隔迫零块对角导频图案，它在每个导频时隙后插入了N_g-1个填零的时隙，这些零所持续的时间大于延迟最长的路径拍所对应的时隙长度，以克服符号间的干扰并使信道的多拍响应转化成多个单拍响应的和。

在通信接收端实施本发明的IRS辅助级联MIMO双稀疏信道状态信息估计方法时，当N_R＝N_I＝N_T＝4和L₁＝2及L₂＝3时，具体实施步骤如下：

步骤一：用角度随机采样对角度域的平面阵列构造1倍阵列尺寸的4×4维字典矩阵D_1R和D_1T以及D_2R和D_2T，其中β₁＝4和β₂＝4；

步骤二：根据多普勒频移确定ψ后得到角度增益向量h的状态空间模型；

步骤三：初始化h(0)为256×1维复随机向量，初始化其估计误差方差为P(0)＝100I₂₅₆，其中I₂₅₆表示256×256维的单位矩阵，设置更新噪声v(s)的协方差矩阵分别为Q＝I₂₅₆；根据信噪比(SNR)确定观测噪声方差其中I₁₆表示16×16维的单位矩阵；

步骤四：通过带有保护间隔迫零块对角化导频，得到4×4维观测矩阵Y(s)并向量化为16×1维向量y(s)；

步骤五：由常规卡尔曼滤波方法计算h(s|s-1)、P(s|s-1)、K(s)、h(s)和P(s)；

步骤六：构造l₁范数的伪测量方程为0＝||h||₁-η(τ)∈，在自适应因子η(τ)的参数设置中，设置η(0)＝0.9，μ＝0.02，T＝10，设置∈的方差r_∈＝40000，将l₁范数的测量方差设置为r＝η²(τ)r_∈，设置h_l(0)和P_l(0)；

步骤七：在当前估计值上计算l₁范数的次梯度，得到1×256的次梯度向量c，然后进行l₁范数修正，依次计算h_l(τ|τ-1)、P_l(τ|τ-1)、k_l(τ)、h_l(τ)和P_l(τ)；

步骤八：判断相邻两次的l₁范数差值的绝对值是否小于10^-3，若小于则停止迭代，若不小于再判断是否到达最大迭代次数20，若达到停止迭代；停止迭代后，将已进行l₁范数迭代修正的h_l(τ)作为当前h(s)，将P_l(τ)作为当前的P(s)，然后进入步骤四至步骤八的由下一个观测向量的递推更新过程；

步骤九：滤波达到稳态后，得到稀疏的角度增益向量h(s)中的角度支集和相应复角度增益，由角度字典得到瞬时的16×1维vec(H(s))和4×4维级联信道H(s)。

在实施例中采用了毫米波级联信道S₁，在类似的实施例：N_R＝N_I＝N_T＝4和L₁＝1及L₂＝1中，采用了毫米波级联信道S₂，在它们的角度字典中所用的方位角和仰角以及相应的角度增益分别如表1和表2所示。依据式(1)、(2)和(3)可生成瞬时的4×4维MIMO级联信道。

表1N_I＝4、L₁＝2、L₂＝3的4×4维MIMO级联信道

表1N_I＝4、L₁＝1、L₂＝1的4×4维MIMO级联信道

在动态场景下，实施例中两个场景S₁和S₂的角度参数不变，真实角度增益的初始值仍然如表1和表2所示，而动态的角度增益按状态空间模型随符号更新。在实施例中，将4个IRS元件在4个连续符号中的相移设置为{1，1，1，1}，{1，j，-1，j}，{1，-1，1，-1}和{1，-j，-1，j}，并随符号s的增加重复这4组相移。

如附图5所示为本发明方法(记为SKF)与未带有自适应因子的稀疏滤波方法(记为GKF)、最优最小二乘(记为Oracle LS)和OMP在静态场景S₁信道参数设置下的归一化均方误差(NMSE)的对比曲线，NMSE定义为：

式中，表示在最后一个导频符号处角度增益向量的估计值，h(S)表示真实角度增益向量值，L表示仿真次数。可以看出，本发明所提出的方法具有最好的NMSE性能，OMP和Oracle LS已失去了对角度支集的估计能力，性能很差；GKF的NMSE也较低。

如附图6所示为本发明方法对其中一个非零值的角度增益h₁的估计跟踪曲线，可以看出，本发明SKF方法在估计过程和估计精度上明显好于一般的GKF方法。

如附图7所示为本发明方法对6个非零角度增益的角度支集和角度增益幅值的估计结果，相较于GKF方法，本发明方法的误差更低。

如附图8所示为在静态场景S₁和SNR＝10dB时，NMSE随导频音观测向量个数的变化曲线，可以看出，在相同的导频音观测向量个数时，本发明方法能获得更好的NMSE性能。

如附图9所示为SNR＝10dB和导频音观测向量数为20时，NMSE随ψ的变化曲线，可以看出，较小ψ时的NMSE较大，随着ψ的增大直到1，NMSE曲线下降，本发明方法SKF在动态场景S₁具有更好的NMSE性能。

如附图10所示为在SNB＝10dB、导频音观测向量数为20和ψ＝0.9时对一个非零的角度增益h₁的动态跟踪误差曲线；在仿真中信道的初始值为表1(a)的数值。可以看出，本发明方法具有较好的动态跟踪效果，能在很大程度上跟踪动态稀疏信道。

如附图11所示为在SNR＝10dB和ψ＝0.9时所估计的角度增益向量的l₁范数随导频音观测向量数的变化曲线，可以看出，本发明方法使用较少的导频音观测向量就可以使l₁范数达到稳态，并且能获得更接近真实l₁范数的估计结果。

如附图12至18所示为信道参数为S₂时的性能曲线，用对附图5-11类似的分析，可得出本发明方法具有优良的IRS辅助级联MIMO双稀疏信道状态信息的估计性能。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (10)

1.一种IRS辅助级联MIMO双稀疏信道状态信息估计方法，其特征在于，方法包括以下步骤：

S1、对于配置N_T根发送天线的通信发送端、带N_I个反射元件的IRS和配置N_R根接收天线的通信接收端所组成的级联传输系统，为通信发送端与IRS端之间的信道H₁，构造接收角度域字典矩阵D_1R和发送角度域字典矩阵D_1T，同时为通信接收端与IRS端之间的信道H₂，构造接收角度域字典矩阵D_2R和发送角度域字典矩阵D_2T，基于上述字典矩阵，得到经角度过采样后的级联信道的表达式，再将该表达式向量化，向量化的表达式中包括待求的级联角度增益向量h；

S2、对S1的级联角度增益向量建立任意的第n个符号对应的状态空间模型h(n)；

S3、初始化级联角度增益向量h为h(0)，初始化级联角度增益向量h的估计误差方差为P(0)，设置状态空间模型中更新噪声的协方差矩阵Q，根据信噪比设置观测噪声方差R；

S4、通过第s个符号的带有保护间隔迫零块对角导频X(s)获得导频音观测矩阵Y(s)，并进行向量化，得到当前观测向量y(s)，所述观测向量y(s)和符号序号s一一对应；

S5、对于当前观测向量，由卡尔曼滤波方法计算角度增益向量的一步预测h(s|s-1)，根据协方差矩阵Q计算一步预测估计方差P(s|s-1)，根据观测噪声方差R计算滤波增益K(s)，根据滤波增益和角度增益向量的一步预测得到当前估计值h(s)，根据滤波增益和一步预测估计方差得到当前估计误差方差P(s)；

S5中第一次采用卡尔曼滤波方法计算时的初始值采用S3的初始化的h(0)和P(0)，其他次数的初始值采用S8的当前第s个符号的级联角度增益向量和估计方差；

S6、构造l₁范数的伪测量方程，设置l₁范数的误差方差，把当前估计值h(s)作为l₁范数迭代修正的初始值h_l(0)，把当前的估计误差方差P(s)作为迭代初始估计误差方差P_l(0)；

S7、用卡尔曼滤波方法计算一步预测h_l(τ|τ-1)，计算一步预测估计误差方差P_l(τ|τ-1)，然后计算迭代滤波增益k_l(τ)，计算当前估计值h_l(τ)和当前估计误差方差P_l(τ)，其中τ为迭代次数；

S7中第一次采用卡尔曼滤波方法计算时的初始值采用S6的初始值h_l(0)和初始估计误差方差P_l(0)；

S8、判断两次相邻的估计值h_l(τ)和h_l(τ-1)的l₁范数的差值绝对值是否小于阈值，若是则停止迭代，反之则判断是否达到最大迭代次数，若是则停止迭代，反之继续执行S7进行迭代；

S8中，停止迭代后，把已进行l₁范数迭代修正的h_l(τ)作为当前第s个符号的级联角度增益向量，把停止迭代时的P_l(τ)作为当前第s个符号的估计误差方差；

S9、更新符号为当前符号s的下一个符号，重复S4～S8，滤波达到稳态后，得到稀疏的角度增益向量中的角度支集和相应的角度增益，根据角度域字典得到瞬时的矩阵向量化的信道，通过矩阵化得到瞬时的级联信道H(s)。

2.根据权利要求1所述的一种IRS辅助级联MIMO双稀疏信道状态信息估计方法，其特征在于，S1的向量化的表达式为：

其中，H表示经角度过采样后的级联信道，g(s)表示IRS在第s个符号上的相移向量，A_R，I和A_I，T分别为在字典导向矢量方向上通信接收端和IRS端之间的角度增益矩阵以及IRS端和通信发送端之间的角度增益矩阵，符号和⊙分别表示克罗内克积和Khatri-Rao积，符号(·)^*，(·)^T和(·)^H分别表示共轭、转置和共轭转置，为待求的级联角度增益向量。

3.根据权利要求2所述的一种IRS辅助级联MIMO双稀疏信道状态信息估计方法，其特征在于，S5的角度增益向量的一步预测h(s|s-1)和一步预测估计方差P(s|s-1)为：

h(s|s-1)＝ψh(s-1)

P(s|s-1)＝ψP(s-1)+(1-ψ²)Q

其中，ψ为状态空间模型中由最大多普勒频移f_D确定的参数，P(s-1)为上一个符号对应的估计误差方差，Q为更新噪声的协方差矩阵。

4.根据权利要求2所述的一种IRS辅助级联MIMO双稀疏信道状态信息估计方法，其特征在于，滤波增益K(s)为：

K(s)＝P(s|s-1)A^H(s)(A(s)P(s|s-1)A^H(s)+R)^-1

其中，R为观测噪声方差，X(s)为带有保护间隔迫零块对角导频符号，(·)^*，(·)^T和(·)^H分别表示共轭、转置和共轭转置；

当前估计值h(s)和当前估计误差方差P(s)为：

h(s)＝h(s|s-1)+K(s)(y(s)-A(s)h(s|s-1))

P(s)＝P(s|s-1)-K(s)A(s)P(s|s-1)

其中，h(s|s-1)为角度增益向量的一步预测，P(s|s-1)为一步预测估计误差方差。

5.根据权利要求1所述的一种IRS辅助级联MIMO双稀疏信道状态信息估计方法，其特征在于，S6的伪测量方程为：

0＝||h(s)||₁-η(τ)∈

其中，η为自适应因子，∈为一个方差固定的随机数，τ为迭代次数。

6.根据权利要求5所述的一种IRS辅助级联MIMO双稀疏信道状态信息估计方法，其特征在于，自适应因子为：

η(τ)＝1-η(0)exp(-μτ/T)

其中，η(0)和μ是形状参数，T为一个固定的正常数，τ为迭代次数。

7.根据权利要求1所述的一种IRS辅助级联MIMO双稀疏信道状态信息估计方法，其特征在于，S7的一步预测h_l(τ|τ-1)和一步预测估计误差方差P_l(τ|τ-1)为：

h_l(τ|τ-1)＝h_l(τ-1)

P_l(τ|τ-1)＝P_l(τ-1)

其中，h_l(τ-1)为上一迭代次数时的估计值，P_l(τ-1)为上一迭代次数时的估计误差方差；

迭代滤波增益k_l(τ)为：

k_l(τ)＝P_l(τ|τ-1)c^H(τ)(c(τ)P_l(τ|τ-1)c^H(τ)+r(τ))^-1

其中，c(τ)为由次梯度向量构成的雅可比矩阵，r(τ)＝η²(τ)r_∈为l₁范数的测量方差，r_∈为随机数∈的固定方差，η为自适应因子。

8.根据权利要求7所述的一种IRS辅助级联MIMO双稀疏信道状态信息估计方法，其特征在于，当前估计值h_l(τ)和当前估计误差方差P_l(τ)为：

h_l(τ)＝h_l(τ|τ-1)+k_l(τ)(0-c(τ)h_l(τ|τ-1))

P_l(τ)＝P_l(τ|τ-1)-k_l(τ)c(τ)P_l(τ|τ-1)

其中，k_l(τ)为迭代滤波增益。

9.根据权利要求1所述的一种IRS辅助级联MIMO双稀疏信道状态信息估计方法，其特征在于，滤波达到稳态的判断条件为：

基于相邻两个符号的S8得到的估计误差方差的主对角线元素和的绝对差是否小于阈值进行判断，若是，则达到稳态，反之未达到稳态。

10.根据权利要求1所述的一种IRS辅助级联MIMO双稀疏信道状态信息估计方法，其特征在于，接收角度域字典矩阵D_1R和发送角度域字典矩阵D_1T分别为：

接收角度域字典矩阵D_2R和发送角度域字典矩阵D_2T分别为：

其中，a_1r，j和a_1t，j，j＝0，...，β₁-1是信道H₁在角度网格点上的导向矢量，a_2r，j和a_2t，j，j＝0，...，β₂-1是信道H₂在角度网格点上的导向矢量。