Estimation of the Pseudo Noise Sequence for Multi-Antenna NPLC-DS-CDMA Signals Based on Tensor Decomposition and Polynomial Library Search

1.   引言

直接序列码分多址(Direct Sequence-Code Division Multiple Access, DS-CDMA)信号因其抗干扰能力强，保密性能好，在军事和民用通信领域中，得到了广泛应用^[1]。在非合作通信系统中，由于接收方接收的是未知伪码序列调制的信号，因此对伪码序列的估计是非合作直扩通信系统中信息获取、解扩及分析的前提和关键。

目前针对DS-CDMA信号伪码序列估计的研究主要集中在短码(Short-Code, SC)^[2]和周期长码(Periodic Long-Code, PLC)^[3,4]调制的DS-CDMA信号。对于非周期长码(Non-Periodic Long-Code, NPLC)DS-CDMA信号，文献[5]将NPLC-DS-CDMA信号建模为含有缺失数据的SC-DS-CDMA信号，然后对缺失数据进行插补并利用嵌套迭代最小二乘投影算法对伪码序列进行估计，但当缺失数据比例较大时，估计性能较差。文献[6]在文献[5]的基础上，对最小二乘法进行了改进，应用在了异步NPLC-DS-CDMA信号伪码序列估计中。随着研究的深入，以张量形式表示的多维数据模型开始应用在DS-CDMA信号的研究中^[7]，文献[8]首先将多天线LC-DS-CDMA信号构建成含有缺失数据的SC-DS-CDMA信号张量模型，然后使用交替三线性分解算法对缺失数据进行插补来估计伪码序列，类似于文献[5]，当缺失数据比例较大时，估计性能较差。文献[9]采用分段的思想，将同步多天线LC-DS-CDMA信号分段构建为SC-DS-CDMA子张量模型，然后利用梯度下降法对各个子张量进行Tucker分解得到伪码序列，估计效果优于插补法。然而，梯度下降法在达到最小点附近的时候收敛速度变慢甚至出现发散情况，并且对初始点的选择极为敏感^[10]，影响估计效果。

本文假设在已经估计出伪码周期，码片速率和用户数的基础上^[11,12]，利用文献[9]的分段思想，将多天线NPLC-DS-CDMA信号建模为SC-DS-CDMA信号子张量，并利用交替最小二乘投影(Alternating Least Squares with Projection, ALSP)算法对每个子张量进行CP(Canonical Polyadic)分解得到唯一的伪码片段因子矩阵和接收增益因子矩阵，然后根据增益矩阵相同的特点估计出置换矩阵，利用置换矩阵和旁瓣能量检测去除伪码片段的排序模糊和相位模糊，得到每个用户的伪码序列。最后，利用多项式库搜索^[13]容错性能好的优点，识别出每个用户伪码序列生成多项式，进一步提高伪码序列估计的正确率。

2.   张量CP分解的信号模型构建

本文采用的基带同步NPLC-DS-CDMA信号多天线接收模型如图1所示，$U$表示用户数，$K$为天线数，${A_u}$, ${b_u}\left( t \right)$, ${s_u}\left( t \right)$分别表示第$u$个用户的信号幅度、信息波形、伪码波形，${a_{ku}}$表示第$k$个天线对第$u$个用户的接收增益，${v_k}\left( t \right)$表示第$k$个接收天线中均值为0、方差为${\sigma ^2}$的高斯白噪声。

图  1  同步NPLC-DS-CDMA信号多天线接收模型

下载: 全尺寸图片 幻灯片

第$k$个天线接收的持续时间为$J$个伪码周期的信号表达式为

其中，$\left\{ {{b_u}\left( m \right) = \pm 1} \right\}_{m = 0}^{M - 1}$为第$u$个用户的信息序列，$q\left( t \right)$为持续时间为${T_{\rm{b}}}$的矩形脉冲，T_b和T分别为信息符号周期和伪码周期，$M$为$J$个伪码周期内所对应的信息符号数目，第$u$个用户一周期的伪码波形${p_u}\left( t \right)$可表示为

其中，$\left\{ {{c_u}\left( l \right) = \pm 1} \right\}_{l = 0}^{L - 1}$为第$u$个用户的伪码序列；$g\left( t \right)$ 表示持续时间为${T_{\rm{c}}}$且具有单位能量的码片波形；${T_{\rm{c}}}$为伪码码片周期；$L = {T/{{T_{\rm{c}}}}}$为伪码序列位数；$G = {{{T_{\rm{b}}}}/{{T_{\rm{c}}}}}$为扩频增益。不失一般性，假设信号以${A_u} = 1$进行等功率传输，在接收端先使信号${y_k}\left( t \right)$通过匹配滤波器$h\left( t \right) = g\left( {{T_{\rm{c}}} - t} \right)$再用码片速率进行采样，则接收信号的离散形式为

其中，$N = JL$表示第$k$个天线接收信号的长度，${y_k}\left( n \right)$的矩阵形式可写为

将矩阵${{{Y}}_k}\left( {k = 1,2, ··· ,K} \right)$按照图2(a)方式排列，可表示为3阶张量${\mathcal{Y}} \in {\mathbb{R}^{L \times J \times K}}$，如图2(b)所示。${\mathcal{Y}}$按照扩频增益$G$可分解为$R = \left\lceil {\dfrac{L}{G}} \right\rceil$个子张量${{\mathcal{Y}}^r} \in {\mathbb{R}^{I \times J \times K}},{\rm{ }}r = 1,2, ··· ,R$，$\left\lceil \cdot \right\rceil$表示向上取整，如图2(c)所示。定义$I$为张量${\mathcal{Y}}$在$L$方向上截取的长度，当伪码周期能被扩频增益整除时，$I = G$；否则

图  2  3阶张量示意图

下载: 全尺寸图片 幻灯片

第$r$个子张量${{\mathcal{Y}}^r} = \left[\kern-0.15em\left[ {y_{ijk}^r} \right]\kern-0.15em\right] \in {\mathbb{R}^{I \times J \times K}}$上的元素可写为

其中，$c_u^r\left( i \right) \!\!=\!\! {c_u}\!\left( {\left( {r \!\! - \!\! 1} \right)G \!\! + \!\! i} \right)$，$b_u^r\left( j \right) \!\!=\!\! {b_u}\!\!\left(\! {\left\lfloor {\left( {j \!\! - \! \! 1} \right)\dfrac{L}{G}} \right\rfloor \!\! +\!\! r} \!\right)$，$v_{ijk}^r$为噪声张量${{{\cal V}}^r} = \left[\kern-0.15em\left[ {v_{ijk}^r} \right]\kern-0.15em\right] \in {\mathbb{R}^{I \times J \times K}}$的元素，$\left\lfloor \cdot \right\rfloor$表示向下取整。则$y_{ijk}^r$表示第$k$个天线接收的第$\left( {\left\lfloor {\left( {j - 1} \right)\dfrac{L}{G}} \right\rfloor + r} \right)$个信息码对应的第$\left( {\left( {r - 1} \right)G + i} \right)$个伪码码片上的数据元素。

若定义因子向量${{c}}_u^r$, ${{b}}_u^r$, ${{a}}_u^r$分别为

则3阶张量${{\mathcal{Y}}^r}$因子分解的元素表达式(6)可以等价表示为

其中，“$\circ$”表示对向量求外积。由式(8)知，3阶张量可分解为$U$个列向量的外积之和，称为张量的CP分解^[14]，如图3所示。

图  3  张量CP分解示意图

下载: 全尺寸图片 幻灯片

图3中，因子向量${{c}}_u^r$, ${{b}}_u^r$, ${{a}}_u^r$$\left( {u = 1,2, ··· ,U} \right)$分别构成张量CP分解的3个因子矩阵，即

其中，因子矩阵${{{C}}^r} \in {\mathbb{R}^{I \times U}}$和${{{B}}^r} \in {\mathbb{R}^{J \times U}}$分别表示$U$个用户第$r$个片段的伪码序列和信息码序列，${{{A}}^r} \in {\mathbb{R}^{K \times U}}$表示$K$个天线对 $U$个用户的接收增益。

3.   伪码序列估计

3.1   基于CP分解ALSP算法的伪码序列片段估计

在实际应用中，为了便于对高维数据进行处理，将张量展开为矩阵形式，根据文献[15]提出的展开方法，对子张量${{\mathcal{Y}}^r}$进行Kolda水平展开为

其中，${{Y}}_{\left( 1 \right)}^r$, ${{Y}}_{\left( 2 \right)}^r$, ${{Y}}_{\left( 3 \right)}^r$分别为模式-1、模式-2、模式-3水平展开矩阵，对应矩阵中的元素分别为$y_{i,\left( {k - 1} \right)J + j}^r$, $y_{j,\left( {k - 1} \right)I + i}^r$, $y_{k,\left( {j - 1} \right)I + j}^r$。${{{Y}}^r}\left( {:,:,k} \right)$和${{{Y}}^r}\left( {:,j,:} \right)$为张量${{\mathcal{Y}}^r}$的正面切片矩阵和侧向切片矩阵。

若不考虑噪声影响，根据平行比例配置剖面原则，3阶张量的正面切片矩阵

则模式-1水平展开矩阵${{Y}}_{\left( 1 \right)}^r$可写作

其中，“$\odot$”表示求两个矩阵的Khatri-Rao积。同理可得，另外两种水平展开矩阵可写作

由式(12)、式(13)和式(14)可知，3阶张量的矩阵化能够由CP分解的3个因子矩阵来表示。由于张量CP分解具有唯一性(相关证明可参考文献[2])，对式(12)、式(13)和式(14)对应的分离优化问题进行求解，得到的因子矩阵分别对应接收信号的伪码片段因子矩阵，信息码片段因子矩阵和接收增益因子矩阵，考虑到噪声的影响，由式(12)可得最小二乘分离优化目标函数为

其中，${\left\| \cdot \right\|_{\rm{F}}}$表示矩阵的Frobenius范数。对应最小二乘解为

其中，“$*$”表示矩阵的Hadamard积，${\left[ \cdot \right]^\dagger }$表示求矩阵的Moore-Penrose逆。同理可得，信息码片段因子矩阵${{{B}}^r}$和接收增益因子矩阵${{{A}}^r}$的最小二乘解为

式(16)中，${{{B}}^r}$中的所有元素${b_{j,u}} \in \left\{ {1, - 1} \right\}$，对式(16)运用最小二乘求解时，首先固定${{{A}}^r}$和${{{B}}^r}$求解${{{C}}^r}$使之满足式(16a)，然后根据${{{A}}^r}$和得到的${{{C}}^r}$求解${{{B}}^r}$使之满足(16b)，最后利用求出的${{{C}}^r}$和${{{B}}^r}$求解${{{A}}^r}$使之满足(16c)，依次循环迭代，直到算法收敛。因此对式(16)的求解实际上采用的是ALSP算法，具体步骤如下：

步骤 1　随机初始化因子矩阵${{C}}_0^r$, ${{B}}_0^r$和${{A}}_0^r$，令迭代次数$\alpha = 0$；

步骤 2　令$\alpha = \alpha + 1$，计算伪码片段因子矩阵${{{C}}^r}$的第$\alpha$次迭代值${{C}}_\alpha ^r$

步骤 3　计算信息码片段因子矩阵${{{B}}^r}$的第$\alpha$次迭代值${{B}}_\alpha ^r$

步骤 4　对${{B}}_\alpha ^r$中的元素进行硬判决，将其映射为$\left\{ { + 1, - 1} \right\}$，即${{B}}_\alpha ^r = {\rm{sign}} \left( {{{B}}_\alpha ^r} \right)$。其中，${\rm{sign}} \left( \cdot \right)$表示符号函数。

步骤 5　计算接收增益因子矩阵${{{A}}^r}$的第$\alpha$次迭代值${{A}}_\alpha ^r$

步骤 6　计算收敛判别函数

步骤 7　判断终止迭代条件，当$\left| \rho \right| < \beta$($\beta$是一个极小的数，通常取${10^{ - 10}}$)时，表示算法收敛，迭代结束；否则，继续执行步骤2至步骤6，直到算法收敛(关于算法收敛性证明可参考文献[5])。

对每个子张量利用CP分解的ALSP算法求解后，可得到$R$个伪码片段因子矩阵，然后对每个伪码片段按照对应用户排序拼接就可得到每个用户的伪码序列。

3.2   基于接收增益矩阵互相关性和旁瓣能量检测的伪码片段组合序列筛选

由3.1节得到的R个伪码片段因子矩阵内部存在排序模糊和相位模糊的问题，由式(7)和式(9)可知，每个子张量的接收增益因子矩阵${{{A}}^r}$相同，${{{A}}^r}\left( {r = 1,2, ··· ,R} \right)$相同位置的列向量表示$K$个天线对同一个用户的接收增益。因此可以利用${{{A}}^r}$对${{{C}}^r}$的列向量重新排序，使得${{{C}}^r}\left( {r = 1,2, ··· ,R} \right)$相同位置的列向量表示同一个用户的伪码片段，完成组合序列的初步筛选。

假设因子矩阵${{{A}}^1}$和${{{A}}^r}$已经估计得到，其互相关矩阵的绝对值为

其中，$\left| \cdot \right|$表示求绝对值。搜索${{{X}}^r}$每行最大值的位置，记为

式(22)表示${{{X}}^r}$第$u$行的最大值位于第${w_u}$列，说明${{{A}}^1}$的第$u$列与${{{A}}^r}$的第${w_u}$列的相关性最大，属于同一个用户，根据这一性质可求出一置换矩阵，记为${{{Z}}^r}$，其中${{{Z}}^r}\left( {u,{w_u}} \right) = 1$，其余位置为0。利用置换矩阵对${{{C}}^r}$进行重新排序，得到排序后的伪码片段矩阵为

解决了排序模糊之后，由于${{{H}}^r}$中每一片段都有正相和反相两种可能，所以每个用户拼接得到的伪码序列(记为${{{H}}_u}$)都有${2^R}$种可能性，通过计算${2^R}$种组合序列的自相关函数的旁瓣能量^[16]，选取旁瓣能量最小的序列即为正确相位的伪码序列。

3.3   基于多项式库搜索的伪码序列估计

由于伪码序列生成多项式个数是有限的，因此可以首先建立一个伪码序列生成多项式库，然后根据生成多项式的阶数划分接收序列为分析矩阵，并构造关于生成多项式系数的方程组，最后在有限的多项式库中进行搜索找到能够与分析矩阵匹配的多项式即为接收序列的生成多项式，进一步提高伪码序列估计的正确率。以m序列为例，根据m序列的线性递推关系${c_i} = \displaystyle \sum\nolimits_{l = 1}^v {{d_l}{c_{i - l}}\left( {i \ge 0} \right)}$，可构造方程组为

其中，${c_i}\left( {i = 0,1, ··· ,T + v} \right) \in {\rm{GF}} \left( 2 \right)$为$v$阶m序列的第$i$位，${d_l}\left( {l = 0,1, ··· ,v} \right) \in {\rm{GF}} \left( 2 \right)$为m序列生成多项式的系数且${d_0} = 1$。符合式(24)方程组的解$\left( {{d_0},{d_1} ··· {d_v}} \right)$即为$v$阶m序列生成多项式$f\left( x \right) = {d_0} + {d_1}{x^1} + ··· + {d_v}{x^v}$。

由于接收序列不可避免地出现误码，则在正确生成多项式下，式(24)中的方程组不一定全部成立，故对$v$阶m序列构造方程组如式(25)。

若接收序列无误码，则${{Z}}$中元素为0的个数为$T + 1$，若接收序列为随机序列，则${{Z}}$中元素为0的个数为${{\left( {T + 1} \right)}/2}$，由于接收序列中含有误码，所以${{Z}}$中元素为0的个数一定小于${{\left( {T + 1} \right)}/2}$，因此取

若$F$越大，则估计的多项式是生成多项式的可能性越大。值得注意的是，若接收序列是m序列的镜像序列(全部反相)，则$\left| F \right|$越大，估计的多项式是生成多项式的可能性越大。

由文献[13]可知，利用多项式库搜索算法对二进制伪随机序列识别时，容错性最大可达到近20%的误码率。因此，对3.2节筛选得到的组合序列进行$- 1 \to 1, \; 1 \to 0$映射后，按照式(25)构造$T + 1$个关于库多项式系数的方程组，($T$可根据实际需要进行设置)，然后根据式(26)计算$F$，得到每个用户伪码序列的生成多项式，从而实现对组合序列纠错，提高伪码序列估计的正确率。

4.   仿真实验与结果分析

不失一般性，实验中信息码采用随机产生的BPSK调制序列，伪码序列采用m序列，${a_{ku}}$是均值为0、方差为1的高斯分布随机数，接收信号长度$N = 1000L$，构造方程组的规模$T = L - v$，信噪比取值范围${\rm{SNR}} = - 20 \sim 0\;{\rm{dB}}$，且每一个信噪比下Monte-Carlo仿真次数为400次。

实验1　算法性能与天线数$K$和用户数$U$的关系。取伪码周期长度$L = 127$，扩频增益$G = 50$，得到用户数$U = 4,6,8,10$在天线数$K = 6,10$的条件下的伪码序列估计性能曲线如图4和图5所示。

图  4  K=6下不同用户数的伪码序列估计性能

下载: 全尺寸图片 幻灯片

图  5  K=10下不同用户数的伪码序列估计性能

下载: 全尺寸图片 幻灯片

由图4和图5可看出，天线数一定情况下，用户数越多，信号越复杂，估计性能越差；对比图4和图5可知，用户数一定的情况下，天线数越多，估计性能越好，这是因为天线数越多，接收的有用信号越多，算法性能越好。且本文算法能够在天线数少于用户数的条件下，信噪比为$- 2\;{\rm{dB}}$时，实现对10个用户伪码序列的无误码估计。

实验2　算法性能与伪码周期长度$L$和扩频增益$G$的关系。取天线数$K = 6$，用户数$U = 4$，得到$L = 127$, $G = 50,\;70$及$L = 511$, $G = 201,\,282$条件下的伪码序列估计性能曲线如图6所示。

图  6  不同L和G下的伪码序列估计性能

下载: 全尺寸图片 幻灯片

由图6可知，在伪码周期长度$L$一定情况下，扩频增益$G$越小，分割的子张量越多，对伪码周期破坏性越大，估计效果越差；在扩频调制比${L/G}$近似相等，即分割的子张量个数相同的情况下，则伪码周期长度$L$越大，伪码序列越接近白噪声，抗干扰性能越好，估计性能越好。

实验3　本文算法同文献[5]及文献[9]算法的伪码序列估计性能对比。取天线数$K = 6$，用户数$U = 4$，伪码周期长度$L = 127$，扩频增益$G = 50, \; 70$，得到不同文献的伪码序列估计性能曲线如图7所示。

图  7  本文算法同文献[5]和文献[9]算法性能对比

下载: 全尺寸图片 幻灯片

由图7可知，本文算法对伪码序列的估计性能优于文献[5]和文献[9]。文献[5]将NPLC-DS-CDMA信号构建成含有缺失数据的信号模型，并对缺失数据进行插补，破坏了原信号的结构信息，估计效果较差，而文献[9]和本文利用张量建模，能够保持各元素内在联系，充分利用信号空间结构信息，因此估计性能要远好于文献[5]；本文利用多项式库搜索对估计序列进行纠错，进一步降低误码率，其估计性能同文献[9]相比有明显提升。

实验4　多项式库搜索对伪码序列估计性能的提升。取天线数$K = 10$，伪码周期长度$L = 127$，扩频增益$G = 50$，得到在不同用户数$U = 4,8$下使用多项式库搜索和未使用多项式库搜索的伪码序列估计性能曲线如图8所示，其中未使用多项式库搜索是指在3.2节后直接计算筛选后的伪码序列的误码率。

图  8  多项式库搜索对伪码序列估计性能提升

下载: 全尺寸图片 幻灯片

由图8可看出，当用户数$U = 4$和$U = 8$时，使用多项式库搜索达到的信噪比容限相对于未使用多项式库搜索分别提升了约$3\;{\rm{dB}}$和$2\;{\rm{dB}}$。这是因为用户数越少，估计性能越好，多项式库搜索在一定容错性下，对序列的纠错性能越好，提升效果越明显。

实验5　本文算法同文献[9]算法的迭代次数对比。取天线数$K = 10$，用户数$U = 4$，码周期长度$L = 127$，扩频增益$G = 50$，得到本文算法和文献[9]算法收敛时的平均迭代次数随信噪比的变化如图9所示。

图  9  本文算法同文献[9]算法迭代次数对比

下载: 全尺寸图片 幻灯片

由图9可看出，在相同信噪比下，本文使用ALSP算法达到收敛时的迭代次数远小于文献[9]使用的梯度下降法，较文献[9]相比具有更快的收敛速度，且迭代次数受信噪比影响较小。

5.   结束语

针对NPLC-DS-CDMA信号伪码估计问题，本文用多天线对信号进行接收构造3阶张量模型并按照扩频增益分解为多个子张量，用ALSP算法对每个张量进行CP分解估计出因子矩阵，提高了算法的收敛速度。然后利用各子张量接收增益矩阵互相关性及m序列旁瓣能量特性筛选出每个用户的伪码序列。最后利用多项式库搜索方法估计出伪码序列的生成多项式，对伪码序列进行纠错，进一步提高伪码序列估计的正确率。仿真实验表明，本文算法能够在天线数为6，信噪比为$- 2\;{\rm{dB}}$时，实现对10个用户伪码序列的无误码估计，且对伪码序列的估计性能及收敛速度与现有算法相比，都有所提高。