毫米波MIMO-OFDM系统中基于3D-MUSIC算法的信道估计方法

毫米波mimo-ofdm系统中基于3d-music算法的信道估计方法
技术领域
1.本发明属于毫米波通信技术领域，尤其涉及一种毫米波mimo-ofdm系统中基于3d-music算法的信道估计方法。


背景技术：

2.毫米波(mm-wave)通信作为未来无线通信系统的关键技术之一，具有通信容量大、短波长、窄波束、抗干扰和安全保密性好等优势，但也正是由于其自身的特点，与传统的微波系统相比，毫米波系统的传播损耗更大。因此，可以在基站(bs)和移动基站(ms)上采用大规模天线阵列来补偿随频率增加的传输损耗，从而保证毫米波通信具有较大的波束成形增益。一种常用的毫米波系统天线阵列处理方法是采用混合多输入多输出(mimo)系统，该系统由模拟射频波束成形器和数字基带mimo处理器组成，其中模拟射频波束成形器可以通过相控阵或透镜阵来有效地实现，且由于毫米波信道的稀疏性，射频(rf)链路数可以显著地小于天线数，这说明混合模拟/数字结构具有较低的实现成本和复杂度。为了设计具有较高波束成形增益的混合mimo系统，需要通过信道估计来获取准确的信道状态信息(csi)。
3.正交频分复用(ofdm)技术因其对抗无线信道频率选择性衰落以及高频谱效率的优点，被广泛应用于宽带无线通信系统。因此，在毫米波mimo-ofdm系统中研究信道估计问题是十分有必要的。由于毫米波信道的脉冲响应仅由少量的散射路径控制，传统的基于丰富散射的mimo信道模型无法表征毫米波信道的空间稀疏性。但是，毫米波信道可以以参数形式进行建模，即以信道路径的离开角(aod)、到达角(aoa)、传输时延和相应的路径增益等参数进行建模。也就是说，毫米波信道估计问题不再是信道矩阵估计问题，而是信道参数估计问题。


技术实现要素：

4.本发明目的在于提供一种毫米波mimo-ofdm系统中基于3d-music算法的信道估计方法,以获取毫米波宽带信道的信道状态信息的技术问题。
5.为解决上述技术问题，本发明的具体技术方案如下：
6.一种毫米波mimo-ofdm系统中基于3d-music算法的信道估计方法，建立上行链路信号传输的宽带稀疏多路径参数化信道模型，其中每条路径由传输时延，到达角(aoa)，离开角(aod)以及路径增益这四个参数指定；通过三维波束空间多信号分类(music)算法同时估计出信道的aoa、时延和aod这三个参数，并通过鸡群优化(cso)算法来确定空间谱上多个谱峰的所在位置；通过最小二乘(ls)法估计出信道的路径增益参数；
7.所述信道估计方法具体包括以下步骤：
8.s1、建立上行链路信号传输的宽带稀疏多路径参数化信道模型，其中每条路径由四个参数指定：传输时延τ
l
，到达角θ
l
，离开角以及路径增益α
l
(t)；
9.s2、通过三维波束空间多信号分类算法同时估计出信道的到达角传输时延和离开角这三个参数；
10.s3、通过结合三维波束空间多信号分类算法和鸡群优化算法来确定空间谱上多个谱峰的所在位置，以同时估计出信道的多路径参数，即传输时延到达角以及离开角
11.s4、通过最小二乘法来估计信道的路径增益参数。
12.进一步的，所述步骤s1具体包括以下步骤：
13.s11、建立宽带稀疏多路径参数化信道模型：接收端bs配有qa根天线，发射端ms配有ua根天线，设置ofdm子载波总数为m0，并选择前m子载波进行信道估计；接收端bs和发射端ms分别采用半波长间距的一维均匀线性阵列，并通过移相器连接到模拟波束成形器和模拟波束成形器连接到qr和ur根射频链，其中qr＜qa，ur＜ua；采用射线跟踪信道模型，即接收端bs和发射端ms之间有l条可解析的物理路径，其中第l条路径由传输时延τ
l
、到达角θ
l
、离开角和路径增益α
l
(t)这四个参数来指定，其中l＝1,
…
,l；划分一个传输时间帧为n个传输时间块，且复增益在不同传输时间块之间是可变的，而到达角θ
l
、离开角和传输时延τ
l
在每个传输帧中是恒定的，即信道变化仅由路径增益引起，则在第n个传输时间块上的时延域时变信道模型为
[0014][0015]
其中n＝1,
…
,n，表示向量的转置，δ(
·
)表示狄垃克函数；
[0016]
在式(1)中对τ进行傅里叶变换，则在第n个传输时间块和第m个子载波上的频域时变信道模型表示为
[0017][0018]
其中m＝1,
…
,m，是采样率为fs的ofdm子载波间距，e是自然常数，是在第m个子载波和第l条路径上的时延扩展，是接收端bs的天线阵列响应，是发射端ms的天线阵列响应，其表达式分别为
[0019]
[0020][0021]
s12、建立上行链路信号传输模型：设置在每个传输时间块中有t0个时隙，其中前t
p
个时隙被用于信道估计，剩余的时隙被用于数据传输；在一个传输时间块期间，设发射端ms的第i条射频链传输的导频信号被定义为满足
[0022][0023]
其中i＝1,
…
,ur，sj是发射端ms的第j条射频链传输的导频信号，其中j＝1,
…
,ur，表示向量sj的共轭转置；设置在每个传输块期间射频链传输的导频信号相同，且所有射频链传输的导频信号对于所有被选定的子载波都是相同的，则所有射频链发送的导频信号为因此在第n个传输时间块和第m个子载波上，接收端波束成形器上的输出信号为
[0024][0025]
其中为噪声矩阵，其每个元素服从f
*
为矩阵f的共轭矩阵，且在接收端bs和发射端ms上的dft波束成形器被分别表示为
[0026][0027][0028]
其中其中
[0029]
根据式(5)中所描述的导频信号的正交性，得到发射端ms的第i条射频链在第n个传输时间块和第m个子载波上的信号，即
[0030][0031]
其中是矩阵f的第i列，是发射端ms的第i条射频链在第n个传输时间块和第m个子载波上的噪声向量，其每个元素服从均值为0方差为σ2的独立高斯分布；将式(2)代入式(9)得到
[0032][0033]
其中
⊙
表示khatri-rao积，
[0034]
是在第m个子载波上的l条路径的时延扩展向量，
[0035]
是l条路径在发射端ms的天线阵列响应矩阵，
[0036]
是l条路径在接收端bs的阵列响应矩阵，
[0037]
是在第n个传输时间块上的l条路径的复增益；
[0038]
将ur条射频链和m个子载波上的向量y
m,i
(n)堆叠在一列中，得到
[0039][0040]
其中
[0041]
是时延扩展矩阵，其第l列向量表示在第l条路径上的m个子载波的时延扩展，被定义为
[0042][0043]
是将ur条射频链和m个子载波上的向量n
m,i
(n)堆叠在一列所得到的向量，其每个元素服从均值为0方差为σ2的独立高斯分布。
[0044]
进一步的，所述步骤s2具体包括以下步骤：
[0045]
s21、输入路径数l；输入信号y(n)；
[0046]
s22、在式(11)中，定义空间矩阵为
[0047][0048]
其第l列表示为
[0049]
s23、计算信号y(n)的协方差矩阵ry，即
[0050][0051]
其中表示计算y(n)y(n)h的期望。
[0052]
s24、对协方差矩阵ry进行特征分解，即
[0053][0054]
其中是以特征值为对角线元素的对角矩阵，矩阵u的每一列是对应于是每一个特征值的特征向量，矩阵是信号子空间，矩阵是噪声子空间，且|bhun|＝0；
[0055]
s25、计算3d-music算法的空间谱
[0056][0057]
其中，τ∈[0,τ
max
]；这里，和分别是模拟波束成形器w和f的波束扇区边界，分别在式(7)和式(8)中被定义；τ
max
是最大时延。
[0058]
s26、搜索3d-music算法空间谱上的l个谱峰，第l个谱峰所对应的坐标即为所估计出的接收端bs与发射端ms之间第l条路径的时延离开角和到达角
[0059]
进一步的，所述步骤s3具体包括以下步骤：
[0060]
s31、初始化：初始化全局最优解集合搜索范围迭代次数为l-1，并令坐标和的适应度值为+∞；
[0061]
s32、计算全局最优：在搜索范围内执行鸡群优化算法，得到全局最优解
[0062]
s33、更新集合将添加到集合中，并对集合升序排序，则集合中适应度值最小的元素即为空间谱上第1个谱峰所对应的位置坐标；
[0063]
s34、迭代开始：令的适应度值为+∞；
[0064]
s35、更新搜索范围：定义其中i为在集合中的索引值，更新搜索范围为和并分别在每个搜索范围内执行鸡群优化算法，输出新的全局最优解和
[0065]
s36、更新集合将和添加到集合中，并对集合升序排序，则集合中适应度值最小的元素即为空间谱上其他谱峰所对应的位置坐标；
[0066]
s37、迭代终止：若达到最大迭代次数，则停止迭代，得到空间谱上l个谱峰所对应的位置坐标；反之，则返回s34继续迭代。
[0067]
进一步的，所述步骤s32和s35中的鸡群优化算法的具体执行步骤如下：
[0068]
步骤1、初始化：包括迭代次数、种群规模、种群更新间隔代数g、鸡群维度d＝3、公鸡母鸡小鸡在鸡群中所占比例、搜索步长和范围；
[0069]
步骤2、计算适应度：根据式(16)的分母计算每个个体的适应度值，并按照适应度值的大小对鸡群个体排序，以确定鸡群的等级制度；其中公鸡的搜索能力最强，处于统治地位，适应度值最小；母鸡的搜索能力稍差，紧跟在公鸡周围搜索食物，适应度值稍大；小鸡的搜索能力最差，只能跟在母鸡周围搜索食物，适应度值最大，可以实现局部搜索功能；
[0070]
步骤3、初始化鸡群：按公鸡个数将鸡群分组，并将母鸡随机分配到各个分组中，以确定公鸡和母鸡的伙伴关系。然后，小鸡随机跟随母鸡，以确定母鸡和小鸡的母子关系。
[0071]
步骤4、迭代开始：若迭代次数是种群更新间隔代数g的倍数，则更新鸡群分组和鸡群中的关系；否则，公鸡、母鸡和小鸡按照各自的位置更新策略来更新自己的位置，并计算更新位置后的适应度值。
[0072]
步骤5、更新全局最优：将新位置的适应度值与原位置的适应度值相比较，如果新位置的适应度值小，就更新个体位置，否则就保持原来的位置不变。
[0073]
步骤6、迭代终止：若达到最大迭代次数，则停止迭代，输出全局最优解(即空间谱上谱峰所对应的位置坐标)；反之，则返回步骤4继续循环迭代。
[0074]
进一步的，所述步骤s4具体包括以下步骤：
[0075]
s41、根据式(13)和估计出的传输时延到达角和离开角得到相应的估计空间矩阵
[0076]
s42、通过最小二乘法估计出接收端bs与发射端ms之间的路径增益，即
[0077][0078]
其中，表示矩阵的逆矩阵。
[0079]
s43、根据式(2)和估计出的传输时延到达角离开角和路径增益估计出接收端bs与发射端ms之间在第n个传输时间块和第m个子载波上的道矩阵，即
[0080][0081]
本发明的毫米波mimo-ofdm系统中基于3d-music算法的信道估计方法，具有以下优点：通过三维波束空间多信号分类(music)算法同时估计出信道的aoa、时延和aod这三个参数，并通过鸡群优化(cso)算法来确定空间谱上多个谱峰的所在位置；通过最小二乘(ls)法估计出信道的路径增益参数。本发明的信道估计方法具有非常高的准确性，以及较好的信道参数估计性能。
附图说明
[0082]
图1为本发明的具体实施方式的方法流程图；
[0083]
图2为本发明的具体实施方式的毫米波混合mimo-ofdm系统示意图；
[0084]
图3(a)为本发明具体实施方式在不同snr下到达角(aoa)的归一化均方误差(nmse)曲线图；
[0085]
图3(b)为本发明具体实施方式在不同snr下离开角(aod)的归一化均方误差(nmse)曲线图；
[0086]
图3(c)为本发明具体实施方式在不同snr下传输时延的归一化均方误差(nmse)曲线图；
[0087]
图3(d)为本发明具体实施方式在不同snr下路径增益的归一化均方误差(nmse)曲线图；
[0088]
图3(e)为本发明具体实施方式在不同snr下信道矩阵的归一化均方误差(nmse)曲线图。
具体实施方式
[0089]
为了更好地了解本发明的目的、结构及功能，下面结合附图，对本发明一种毫米波mimo-ofdm系统中基于3d-music算法的信道估计方法做进一步详细的描述。
[0090]
如图2所示，考虑毫米波混合mimo-ofdm系统，接收端bs配有qa根天线，发射端ms配有ua根天线。假设ofdm子载波总数为m0，并选择前m子载波进行信道估计。bs和ms分别采用半波长间距的一维均匀线性阵列(ula)，并通过移相器连接到模拟波束成形器和模拟波束成形器连接到qr(＜qa)和ur(＜ua)根射频(rf)链。
[0091]
如图1所示，本发明包括以下步骤：
[0092]
s1.建立上行链路信号传输的宽带稀疏多路径参数化信道模型，其中每条路径由四个参数指定：传输时延τ
l
，到达角(aoa)θ
l
，离开角以及路径增益α
l
(t)。所述步骤s1具体包括：
[0093]
s11.建立宽带稀疏多路径参数化信道模型：考虑毫米波信道的稀疏性，本文采用射线跟踪信道模型，即bs和ms之间有l条可解析的物理路径，其中第l(l＝1,
…
,l)条路径由传输时延τ
l
、到达角(aoa)θ
l
、离开角和路径增益α
l
(t)这四个参数来指定。划分一个传输时间帧为n个传输时间块，且复增益在不同传输块之间是可变的，而aoa、aod和传输延迟在每个传输帧中是恒定的，即信道变化仅由路径增益引起，则在第n(n＝1,
…
,n)个传输时间块上的时延域时变信道模型为
[0094][0095]
其中表示向量的转置，δ(
·
)表示狄垃克函数；
[0096]
在式(1)中对τ进行傅里叶变换，则在第n个传输时间块和第m(m＝1,
…
,m)个子载波上的频域时变信道模型表示为
[0097][0098]
其中是采样率为fs的ofdm子载波间距，e是自然常数，是在第m个子载波和第l条路径上的时延扩展，和分别是bs和ms的天线阵列响应，其表达式分别为
[0099][0100][0101]
s12.建立上行链路信号传输模型：假设在每个传输时间块中有t0个时隙，其中前t
p
个时隙被用于信道估计，剩余的时隙被用于数据传输。在一个传输时间块期间，设发射端ms的第i(i＝1,
…
,ur)条rf链传输的导频信号被定义为满足
[0102][0103]
其中sj是发射端ms的第j(j＝1,
…
,ur)条射频链传输的导频信号，表示向量sj的共轭转置。设置在每个传输块期间rf链传输的导频信号相同，且所有rf链传输的导频信号对于所有被选定的子载波都是相同的，则所有rf链发送的导频信号为因此在第n个传输时间块和第m个子载波上，接收端波束成形器上的输出信号为
[0104][0105]
其中为噪声矩阵，其每个元素服从f
*
为矩阵f的共轭矩阵，且在接收端bs和发射端ms上的dft波束成形器被分别表示为
[0106][0107][0108]
其中其中
[0109]
根据式(5)中所描述的导频信号的正交性，可以得到ms的第i条rf链在第n个传输时间块和第m个子载波上的信号，即
[0110][0111]
其中是矩阵f的第i列，是发射端ms的第i条射频链在第n个传输时间块和第m个子载波上的噪声向量，其每个元素服从均值为0方差为σ2的独立高斯分布；将式(2)代入式(9)可以得到
[0112][0113]
其中
⊙
表示khatri-rao积，
[0114]
是在第m个子载波上的l条路径的时延扩展向量，
[0115]
是l条路径在发射端ms的天线阵列响应矩阵，
[0116]
是l条路径在接收端bs的阵列响应矩阵，
[0117]
是在第n个传输时间块上的l条路径的复增益。
[0118]
将ur条射频链和m个子载波上的向量y
m,i
(n)堆叠在一列中，得到
[0119][0120]
其中
[0121]
是时延扩展矩阵，其第l列向量表示在第l条路径上的m个子载波的时延扩展，被定义为
[0122][0123]
是将ur条射频链和m个子载波上的向量n
m,i
(n)堆叠在一列所得到的向量，其每个元素服从均值为0方差为σ2的独立高斯分布。
[0124]
s2.通过三维波束空间music算法同时估计出信道的aoa、时延和aod这三个参数。所述步骤s2具体包括：
[0125]
s21：输入路径数l；输入信号y(n)。
[0126]
s22：在式(11)中，定义空间矩阵为
[0127][0128]
其第l列可以表示为
[0129]
s23：计算信号y(n)的协方差矩阵ry，即
[0130][0131]
其中表示计算y(n)y(n)h的期望。
[0132]
s24：对协方差矩阵ry进行特征分解，即
[0133]
[0134]
其中是以特征值为对角线元素的对角矩阵，矩阵u的每一列是对应于是每一个特征值的特征向量，矩阵us是信号子空间，矩阵是噪声子空间，且|bhun|＝0；
[0135]
s25：计算3d-music算法的空间谱
[0136][0137]
其中，τ∈[0,τ
max
]。这里，和分别是模拟波束成形器w和f的波束扇区边界，分别在式(7)和式(8)中被定义；τ
max
是最大时延。
[0138]
s26：搜索3d-music算法空间谱上的l个谱峰，第l个谱峰所对应的坐标即为所估计出的bs与ms之间第l条路径的时延和
[0139]
s3.通过结合三维波束空间music算法和cso算法来确定空间谱上多个谱峰的所在位置，以同时估计出信道的多路径参数。所述步骤s3具体包括：
[0140]
s31.初始化：初始化全局最优解集合搜索范围迭代次数为l-1，并令坐标和的适应度值为+∞；
[0141]
s32.计算全局最优：在搜索范围内执行cso算法，得到全局最优解
[0142]
s33.更新集合将添加到集合中，并对集合升序排序，则集合中适应度值最小的元素即为空间谱上第1个谱峰所对应的位置坐标。
[0143]
s34.迭代开始：令的适应度值为+∞。
[0144]
s35.更新搜索范围：定义其中i为在集合中的索引值，更新搜索范围为和并分别在每个搜索范围内执行cso算法，输出新的全局最优解和
[0145]
s36.更新集合将和添加到集合中，并对集合升序排序，则集合中适应度值最小的元素即为空间谱上其他谱峰所对应的位置坐标。
[0146]
s37.迭代终止：若达到最大迭代次数，则停止迭代，得到空间谱上l个谱峰所对应的位置坐标；反之，则返回s34继续迭代。
[0147]
所述步骤s32和s35中的鸡群优化算法的具体执行步骤如下：
[0148]
步骤1、初始化：初始化迭代次数、种群规模、种群更新间隔代数g、鸡群维度d＝3、
公鸡母鸡小鸡在鸡群中所占比例、搜索步长和范围等。
[0149]
步骤2、计算适应度：根据式(16)的分母计算每个个体的适应度值，并按照适应度值的大小对鸡群个体排序，以确定鸡群的等级制度。其中公鸡的搜索能力最强，处于统治地位，适应度值最小；母鸡的搜索能力稍差，紧跟在公鸡周围搜索食物，适应度值稍大；小鸡的搜索能力最差，只能跟在母鸡周围搜索食物，适应度值最大，可以实现局部搜索功能。
[0150]
步骤3、初始化鸡群：按公鸡个数将鸡群分组，并将母鸡随机分配到各个分组中，以确定公鸡和母鸡的伙伴关系。然后，小鸡随机跟随母鸡，以确定母鸡和小鸡的母子关系。
[0151]
步骤4、迭代开始：若迭代次数是种群更新间隔代数g的倍数，则更新鸡群分组和鸡群中的关系；否则，公鸡、母鸡和小鸡按照各自的位置更新策略来更新自己的位置，并计算更新位置后的适应度值。
[0152]
步骤5、更新全局最优：将新位置的适应度值与原位置的适应度值相比较，如果新位置的适应度值小，就更新个体位置，否则就保持原来的位置不变。
[0153]
步骤6、迭代终止：若达到最大迭代次数，则停止迭代，输出全局最优解(即空间谱上谱峰所对应的位置坐标)；反之，则返回步骤4继续循环迭代。
[0154]
s4.通过最小二乘法来估计信道的路径增益参数，所述步骤s4具体包括：
[0155]
s41.根据式(13)和估计出的传输时延到达角和离开角可以得到相应的估计空间矩阵
[0156]
s42.通过最小二乘(ls)法可以估计出bs与ms之间的路径增益，即
[0157][0158]
其中，表示矩阵的逆矩阵。
[0159]
s43.根据式(2)和估计出的传输时延到达角离开角和路径增益可以估计出bs与ms之间在第n个传输时间块和第m个子载波上的信道矩阵，即
[0160][0161]
图3为本发明所提出的信道估计方法在不同snr下，传输时延、aoa、aod、路径增益以及信道矩阵的归一化均方误差(nmse)曲线图。该小区由一个配有16根天线和8个rf链路的bs，以及多个配有8根天线和4个rf链路的ms组成。设bs与ms之间的信道有2条物理路径，每个传输时间帧被划分为100个传输时间块，且选择前8个子载波进行上行信道估计。可以看出，采用本具体实施方式所估计出的传输时延、aoa、aod、路径增益以及信道矩阵具有非常低的nmse，即本发明所提出的信道估计方法具有非常高的准确性。
[0162]
可以理解，本发明是通过一些实施例进行描述的，本领域技术人员知悉的，在不脱离本发明的精神和范围的情况下，可以对这些特征和实施例进行各种改变或等效替换。另外，在本发明的教导下，可以对这些特征和实施例进行修改以适应具体的情况及材料而不会脱离本发明的精神和范围。因此，本发明不受此处所公开的具体实施例的限制，所有落入本技术的权利要求范围内的实施例都属于本发明所保护的范围内。